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

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

Revision: 4137
          http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4137&view=rev
Author:   jswhit
Date:     2007-11-06 15:03:21 -0800 (Tue, 06 Nov 2007)

Log Message:
-----------
fix Antarctica for 'geos' projection, use npy.sum not python built-in, 
refine checks for bad polygons to avoid 'Bus errors' in ctypes libgeos
interface.

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-06 21:33:37 UTC (rev 4136)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py	2007-11-06 23:03:21 UTC (rev 4137)
@@ -792,9 +792,9 @@
                 # coordinates (this saves time, especially for small map
                 # regions and high-resolution boundary geometries).
                 if not containsPole:
-                    # close Antarctica for cylindrical and 
-                    # pseudo-cylindrical projections.
-                    if name == 'gshhs' and self.projection in ['cyl','merc','mill','moll','robin','sinu']:
+                    # close Antarctica for projections in which this
+                    # is necessary.
+                    if name == 'gshhs' and self.projection in ['cyl','merc','mill','moll','robin','sinu','geos']:
                         b = npy.asarray(poly.boundary)
                         lons = b[:,0]
                         lats = b[:,1]
@@ -818,10 +818,7 @@
                     # if polygon instersects map projection
                     # region, process it.
                     if poly.intersects(self._boundarypolyll):
-                        try:
-                            poly = poly.intersection(self._boundarypolyll)
-                        except:
-                            continue
+                        poly = poly.intersection(self._boundarypolyll)
                         # create iterable object with geometries
                         # that intersect map region.
                         if hasattr(poly,'geoms'): 
@@ -866,10 +863,11 @@
                         bx, by = maptran(blons, blats)
                     else:
                         bx, by = self(blons, blats)
-                    mask = npy.logical_and(bx<1.e20,by<1.e20)
+                    goodmask = npy.logical_and(bx<1.e20,by<1.e20)
+                    badmask = npy.logical_or(bx>1.e20,by>1.e20)
                     # if less than two points are valid in 
                     # map proj coords, skip this geometry.
-                    if sum(mask) <= 1: continue
+                    if npy.sum(goodmask) <= 1: continue
                     if name == 'gshhs':
                         # create a polygon object for coastline
                         # geometry.
@@ -878,8 +876,8 @@
                         # if not a polygon,
                         # just remove parts of geometry that are undefined
                         # in this map projection.
-                        bx = npy.compress(mask, bx)
-                        by = npy.compress(mask, by)
+                        bx = npy.compress(goodmask, bx)
+                        by = npy.compress(goodmask, by)
                         # for orthographic projection, all points
                         # outside map projection region are eliminated
                         # with the above step, so we're done.
@@ -891,11 +889,15 @@
                         poly = LineShape(zip(bx,by))
                     # if polygon instersects map projection
                     # region, process it.
-                    if boundarypolyxy.intersects(poly):
-                        try:
-                            poly = boundarypolyxy.intersection(poly)
-                        except:
-                            continue
+                    if not npy.sum(badmask) and boundarypolyxy.intersects(poly):
+                        #print poly.is_valid, poly.geom_type
+                        #import pylab
+                        #a = npy.asarray(boundarypolyxy.boundary)
+                        #b = npy.asarray(poly.boundary)
+                        #pylab.plot(a[:,0],a[:,1],'b')
+                        #pylab.plot(b[:,0],b[:,1],'b')
+                        #pylab.show()
+                        poly = boundarypolyxy.intersection(poly)
                         # create iterable object with geometries
                         # that intersect map region.
                         if hasattr(poly,'geoms'): 
@@ -1131,10 +1133,10 @@
             test2 = npy.fabs(xa-self.llcrnrx) < delx
             test3 = npy.fabs(ya-self.urcrnry) < dely
             test4 = npy.fabs(ya-self.llcrnry) < dely
-            hasp1 = sum(test1*test3)
-            hasp2 = sum(test2*test3)
-            hasp4 = sum(test2*test4)
-            hasp3 = sum(test1*test4)
+            hasp1 = npy.sum(test1*test3)
+            hasp2 = npy.sum(test2*test3)
+            hasp4 = npy.sum(test2*test4)
+            hasp3 = npy.sum(test1*test4)
             if not hasp1 or not hasp2 or not hasp3 or not hasp4:
                 xy = zip(xa.tolist(),ya.tolist())
                 if self.coastpolygontypes[np] != 2:


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