Hi,

the Basemap allskymap example (examples/allskymap.py in current git and 
basemap rev 1.0.2) plots buggy tissots if they cross the zero meridian and/or 
the poles (see http://web.physik.rwth-
aachen.de/~winchen/unfixed_polarRegions_llskymap.png ).

The attached patch fixes this issue (see http://web.physik.rwth-
aachen.de/~winchen/fixed_polarRegions_llskymap.png ).

Best regards.


Tobias

From ff382103981d27a5470eecf727fe962208a505e6 Mon Sep 17 00:00:00 2001
From: Tobias Winchen <winc...@physik.rwth-aachen.de>
Date: Wed, 25 Jan 2012 09:14:12 +0100
Subject: [PATCH 1/2] Added test cases

---
 examples/allskymap.py |    8 ++++++++
 1 files changed, 8 insertions(+), 0 deletions(-)

diff --git a/examples/allskymap.py b/examples/allskymap.py
index 1dc32b0..b06d70b 100644
--- a/examples/allskymap.py
+++ b/examples/allskymap.py
@@ -377,6 +377,14 @@ if __name__ == '__main__':
     # lower right:
     poly = map.tissot(185, -40, 10, 100)
 
+    #center 
+    poly = map.tissot(0, 0, 10, 100, lw=1, color='m')
+    poly = map.tissot(0, -30, 10, 100, lw=1, color='m')
+
+    # poles
+    poly = map.tissot(120, -83, 10, 100, lw=1, color='m')
+    poly = map.tissot(0, 84, 10, 100, lw=1, color='m')
+
     # Plot the tissot centers as "+" symbols.  Note the top left symbol
     # would cross the limb without the clip_path argument; this might be
     # desired to enhance visibility.
-- 
1.7.8.3


From 1dbabb94b26b157e35b9500fadc3681642433c3a Mon Sep 17 00:00:00 2001
From: Tobias Winchen <winc...@physik.rwth-aachen.de>
Date: Wed, 25 Jan 2012 11:58:26 +0100
Subject: [PATCH 2/2] Fix for polar regions and regions crossing zero meridian

---
 examples/allskymap.py |   58 +++++++++++++++++++++++++++++++++++++++++++++++-
 1 files changed, 56 insertions(+), 2 deletions(-)

diff --git a/examples/allskymap.py b/examples/allskymap.py
index b06d70b..ff93cc7 100644
--- a/examples/allskymap.py
+++ b/examples/allskymap.py
@@ -251,10 +251,15 @@ class AllSkyMap(Basemap):
         returns a list of matplotlib.patches.Polygon objects, with two polygons
         when the tissot crosses the limb, and just one polygon otherwise.
         """
-        
+
+        #only use modified version if we are in a critical region
+        if ((self.east_hem(lon_0) and self.east_hem(lon_0+radius_deg)) or (not self.east_hem(lon_0) and not
+            self.east_hem(lon_0-radius_deg))) and abs(lat_0) + radius_deg <90:
+          return Basemap.tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs)
+
         # TODO:  Just return the polygon (not a list) when there is only one
         # polygon?  Or stick with the list for consistency?
-        
+
         # This is based on Basemap.tissot, but addresses a limitation of that
         # method by handling tissots that cross the limb of the map by finding
         # separate polygons in the eastern and western hemispheres comprising
@@ -268,6 +273,55 @@ class AllSkyMap(Basemap):
         delaz = 360./npts
         az = az12
         last_lon = lon_0
+
+        #handling of the poles
+        if (abs(lat_0) + radius_deg >= 90 ):
+            # Use half of the points for the part inside the map, the
+            # other half of the points is at the map border
+            lats = zeros(npts/2)
+            lons = zeros(npts/2)
+            for n in range(npts/2):
+              az = az+delaz *2
+              lon, lat, az21 = g.fwd(lon_0, lat_0, az, dist)
+              lons[n] = lon
+              lats[n] = lat
+
+            a = list(argsort(lons))
+            lons = lons[a]
+            lats = lats[a]
+            x1,y1 = self(lons, lats)
+
+            # Half of the remaining points are used on eastern and
+            # western hemisphere.
+            N = npts/4
+            segs = []
+            dL = (90-abs(lats[0]))/ (N-1)
+            r = range(N)
+
+            # For the south-pole the order of points is changed to plot
+            # the correct polygon
+            if lat_0 < 0:
+              r.reverse()
+              segs.extend(zip(x1,y1))
+            #First Half of the map border
+            x, y = self(-180 * sign(lat_0) *ones(N) , sign(lat_0)* (90 - array(r)* dL) )
+            segs.extend(zip(x,y))
+
+            if lat_0 > 0:
+              segs.extend(zip(x1,y1))
+
+            #Second half of the map border
+            r.reverse()
+            x, y = self(180 * sign(lat_0) *ones(N) , sign(lat_0)* (90 - array(r)* dL) )
+            segs.extend(zip(x,y))
+
+            #z = array(segs)
+            #self.plot(z[:,0], z[:,1])
+            poly = Polygon(segs, **kwargs)
+            ax.add_patch(poly)
+            return [poly]
+
+
         # Note adjacent and opposite edge longitudes, in case the tissot
         # runs over the edge.
         if start_hem:  # eastern case
-- 
1.7.8.3

Attachment: signature.asc
Description: This is a digitally signed message part.

------------------------------------------------------------------------------
Keep Your Developer Skills Current with LearnDevNow!
The most comprehensive online learning library for Microsoft developers
is just $99.99! Visual Studio, SharePoint, SQL - plus HTML5, CSS3, MVC3,
Metro Style Apps, more. Free future releases when you subscribe now!
http://p.sf.net/sfu/learndevnow-d2d
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

Reply via email to