SF.net SVN: matplotlib: [4201] trunk/toolkits/basemap-testing/lib/ matplotlib/toolkits/basemap/basemap.py

[<js...@us...>]

Revision: 4201
          http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4201&view=rev
Author:   jswhit
Date:     2007-11-09 14:11:14 -0800 (Fri, 09 Nov 2007)

Log Message:
-----------
check to see if Greenwich meridian is in map region.  If not, some
speedups can be had.

Modified Paths:
--------------
    trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py

Modified: trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py
===================================================================
--- trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py	2007-11-09 21:28:58 UTC (rev 4200)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py	2007-11-09 22:11:14 UTC (rev 4201)
@@ -775,6 +775,15 @@
             b = npy.asarray(self._boundarypolyll.boundary)
             blons = b[:,0]; blats = b[:,1]
             boundarypolyxy = PolygonShape(zip(*maptran(blons,blats)))
+        # for map regions that don't contain Pole, see if Greenwich
+        # meridian is contained.
+        if not containsPole:
+            lon_0,lat_0 = self(self.xmin+0.5*(self.xmax-self.xmin),\
+                               self.ymin+0.5*(self.ymax-self.ymin),inverse=True)
+            GM0 = PointShape((0,lat_0))
+            GM360 = PointShape((360,lat_0))
+            hasGM0 = GM0.within(boundarypolyll)
+            hasGM360 = GM360.within(boundarypolyll)
         for line in bdatmetafile:
             linesplit = line.split()
             area = float(linesplit[1])
@@ -830,11 +839,16 @@
                     # (so as to properly treat polygons that cross 
                     # Greenwich meridian).
                     if not antart:
-                        blons = b[:,0]-360
-                        poly1 = Shape(zip(blons,blats))
-                        blons = b[:,0]+360
-                        poly2 = Shape(zip(blons,blats))
-                        polys = [poly1,poly,poly2]
+                        if hasGM0:
+                            blons = b[:,0]-360
+                            poly1 = Shape(zip(blons,blats))
+                            polys = [poly1,poly]
+                        elif hasGM360:
+                            blons = b[:,0]+360
+                            poly2 = Shape(zip(blons,blats))
+                            polys = [poly,poly2]
+                        else:
+                            polys = [poly]
                     else: # Antartica already extends from -360 to +720.
                         polys = [poly]
                     for poly in polys:
@@ -906,16 +920,6 @@
                     # if geometry instersects map projection
                     # region, and doesn't have any invalid points, process it.
                     if not badmask.any() and boundarypolyxy.intersects(poly):
-                        #if not poly.is_valid:
-                        #    print poly.geom_type, poly.is_ring, boundarypolyll.is_valid
-                        #    import pylab
-                        #    a = npy.asarray(boundarypolyxy.boundary)
-                        #    b = npy.asarray(poly.boundary)
-                        #    poly2 = boundarypolyxy.intersection(poly)
-                        #    c = npy.asarray(poly2.boundary)
-                        #    pylab.plot(a[:,0],a[:,1],'b')
-                        #    pylab.fill(c[:,0],c[:,1],'r')
-                        #    pylab.show()
                         # if geometry intersection calculation fails,
                         # just move on.
                         try:


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.