SF.net SVN: matplotlib:[7932] trunk/toolkits/basemap

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

Revision: 7932
          http://matplotlib.svn.sourceforge.net/matplotlib/?rev=7932&view=rev
Author:   jswhit
Date:     2009-11-04 16:07:22 +0000 (Wed, 04 Nov 2009)

Log Message:
-----------
added is_land method to check if a point is over land or water.

Modified Paths:
--------------
    trunk/toolkits/basemap/Changelog
    trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py

Modified: trunk/toolkits/basemap/Changelog
===================================================================
--- trunk/toolkits/basemap/Changelog	2009-11-03 21:17:48 UTC (rev 7931)
+++ trunk/toolkits/basemap/Changelog	2009-11-04 16:07:22 UTC (rev 7932)
@@ -1,4 +1,6 @@
 version 0.99.5 (not yet released)
+	   * add "is_land" method to check whether a point is over land or
+             water.
 	   * geos-3.1.1 now required. 3.1.1 source included (instead of 2.2.3).
            * shiftgrid no longer requires a cyclic point to be present
              (patch from Eric Bruning).

Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2009-11-03 21:17:48 UTC (rev 7931)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2009-11-04 16:07:22 UTC (rev 7932)
@@ -801,6 +801,20 @@
                 else:
                     coastsegs.append(seg)
             self.coastsegs = coastsegs
+        # create geos Polygon structures for land areas.
+        # currently only used in is_land method.
+        self.landpolygons=[]
+        self.lakepolygons=[]
+        #self.islandinlakepolygons=[]
+        #self.lakeinislandinlakepolygons=[]
+        x, y = zip(*self.coastpolygons)
+        for x,y,type in zip(x,y,self.coastpolygontypes):
+            b = np.asarray([x,y]).T
+            if type == 1: self.landpolygons.append(_geoslib.Polygon(b))
+            if type == 2: self.lakepolygons.append(_geoslib.Polygon(b))
+            #if type == 3: self.islandinlakepolygons.append(_geoslib.Polygon(b))
+            #if type == 4: self.lakeinislandinlakepolygons.append(_geoslib.Polygon(b))
+
     # set __init__'s docstring
     __init__.__doc__ = _Basemap_init_doc
 
@@ -1537,6 +1551,23 @@
         self.set_axes_limits(ax=ax)
         return rivers
 
+    def is_land(self,xpt,ypt):
+        """
+        Returns True if the given x,y point (in projection coordinates) is
+        over land, False otherwise.  The definition of land is based upon
+        the GSHHS coastline polygons associated with the class instance.
+        Points over lakes inside land regions are not counted as land points.
+        """
+        landpt = False
+        for poly in self.landpolygons:
+            landpt = _geoslib.Point((xpt,ypt)).within(poly)
+            if landpt: break
+        lakept = False
+        for poly in self.lakepolygons:
+            lakept = _geoslib.Point((xpt,ypt)).within(poly)
+            if lakept: break
+        return landpt and not lakept
+
     def readshapefile(self,shapefile,name,drawbounds=True,zorder=None,
                       linewidth=0.5,color='k',antialiased=1,ax=None):
         """


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