From: <js...@us...> - 2009-11-04 16:07:33
|
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. |