Menu
▾
▴
matplotlib-checkins
|
From: <js...@us...> - 2007-11-01 03:18:13
|
Revision: 4076
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4076&view=rev
Author: jswhit
Date: 2007-10-31 20:18:09 -0700 (Wed, 31 Oct 2007)
Log Message:
-----------
initial Shapely integration
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-01 03:16:07 UTC (rev 4075)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-01 03:18:09 UTC (rev 4076)
@@ -17,6 +17,8 @@
from matplotlib.numerix.mlab import squeeze
from matplotlib.cbook import popd, is_scalar
from shapelib import ShapeFile
+from shapely.geometry import Polygon as PolygonShape
+from shapely import wkb
# basemap data files now installed in lib/matplotlib/toolkits/basemap/data
basemap_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
@@ -676,521 +678,23 @@
area_thresh = 10.
else:
raise ValueError, "boundary resolution must be one of 'c','l','i' or 'h'"
- # read in coastline data (only those polygons whose area > area_thresh).
- coastlons = []; coastlats = []; coastsegind = []; coastsegtype = []
- msg = """
-Unable to open boundary dataset file. Only the 'crude', 'low'
-and 'intermediate' resolution datasets are installed by default. If you
-are requesting a 'high' resolution dataset, you need to download
-and install those files manually (see the basemap README for details)."""
- try:
- bdatfile = open(os.path.join(basemap_datadir,'gshhs_'+resolution+'.txt'))
- except:
- raise IOError, msg
- for line in bdatfile:
- linesplit = line.split()
- if line.startswith('P'):
- area = float(linesplit[5])
- west,east,south,north = float(linesplit[6]),float(linesplit[7]),float(linesplit[8]),float(linesplit[9])
- typ = int(linesplit[3])
- useit = self.latmax>=south and self.latmin<=north and area>area_thresh
- if useit:
- coastsegind.append(len(coastlons))
- coastsegtype.append(typ)
- continue
- # lon/lat
- if useit:
- lon, lat = [float(val) for val in linesplit]
- coastlons.append(lon)
- coastlats.append(lat)
- coastsegtype.append(typ)
- coastsegind.append(len(coastlons))
-
- # read in country boundary data.
- cntrylons = []; cntrylats = []; cntrysegind = []
- try:
- bdatfile = open(os.path.join(basemap_datadir,'countries_'+resolution+'.txt'))
- except:
- raise IOError, msg
- for line in bdatfile:
- linesplit = line.split()
- if line.startswith('>'):
- west,east,south,north = float(linesplit[7]),float(linesplit[8]),float(linesplit[9]),float(linesplit[10])
- useit = self.latmax>=south and self.latmin<=north
- if useit: cntrysegind.append(len(cntrylons))
- continue
- # lon/lat
- if useit:
- lon, lat = [float(val) for val in linesplit]
- cntrylons.append(lon)
- cntrylats.append(lat)
- cntrysegind.append(len(cntrylons))
-
- # read in state boundaries (Americas only).
- statelons = []; statelats = []; statesegind = []
- try:
- bdatfile = open(os.path.join(basemap_datadir,'states_'+resolution+'.txt'))
- except:
- raise IOError, msg
- for line in bdatfile:
- linesplit = line.split()
- if line.startswith('>'):
- west,east,south,north = float(linesplit[7]),float(linesplit[8]),float(linesplit[9]),float(linesplit[10])
- useit = self.latmax>=south and self.latmin<=north
- if useit: statesegind.append(len(statelons))
- continue
- # lon/lat
- if useit:
- lon, lat = [float(val) for val in linesplit]
- statelons.append(lon)
- statelats.append(lat)
- statesegind.append(len(statelons))
-
- # read in major rivers.
- riverlons = []; riverlats = []; riversegind = []
- try:
- bdatfile = open(os.path.join(basemap_datadir,'rivers_'+resolution+'.txt'))
- except:
- raise IOError, msg
- for line in bdatfile:
- linesplit = line.split()
- if line.startswith('>'):
- west,east,south,north = float(linesplit[7]),float(linesplit[8]),float(linesplit[9]),float(linesplit[10])
- useit = self.latmax>=south and self.latmin<=north
- if useit: riversegind.append(len(riverlons))
- continue
- # lon/lat
- if useit:
- lon, lat = [float(val) for val in linesplit]
- riverlons.append(lon)
- riverlats.append(lat)
- riversegind.append(len(riverlons))
-
- # extend longitudes around the earth a second time
- # so valid longitudes can range from -360 to 720.
- # This means a lot of redundant processing is done when
- # creating the class instance, but it a lot easier to figure
- # out what to do when the projection domain straddles the
- # Greenwich meridian.
- coastlons2 = [lon+360. for lon in coastlons]
- cntrylons2 = [lon+360. for lon in cntrylons]
- statelons2 = [lon+360. for lon in statelons]
- riverlons2 = [lon+360. for lon in riverlons]
- coastlons3 = [lon-360. for lon in coastlons]
- cntrylons3 = [lon-360. for lon in cntrylons]
- statelons3 = [lon-360. for lon in statelons]
- riverlons3 = [lon-360. for lon in riverlons]
-
- # transform coastline polygons to native map coordinates.
- xc,yc = proj(NX.array(coastlons),NX.array(coastlats))
- xc = xc.tolist(); yc = yc.tolist()
- xc2,yc2 = proj(NX.array(coastlons2),NX.array(coastlats))
- xc3,yc3 = proj(NX.array(coastlons3),NX.array(coastlats))
- xc2 = xc2.tolist(); yc2 = yc2.tolist()
- xc3 = xc3.tolist(); yc3 = yc3.tolist()
-
- # set up segments in form needed for LineCollection,
- # ignoring 'inf' values that are off the map.
- segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:])]
- segmentsll = [zip(coastlons[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:])]
- segtypes = [i for i in coastsegtype[:-1]]
- segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
- segmentsll2 = [zip(coastlons2[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
- segtypes2 = [i for i0,i1,i in zip(coastsegind[:-1],coastsegind[1:],coastsegtype[:-1]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
- segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
- segmentsll3 = [zip(coastlons3[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
- segtypes3 = [i for i0,i1,i in zip(coastsegind[:-1],coastsegind[1:],coastsegtype[:-1]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
- self.coastsegs = segments+segments2+segments3
- self.coastsegsll = segmentsll+segmentsll2+segmentsll3
- self.coastsegtypes = segtypes+segtypes2+segtypes3
-
- # same as above for country segments.
- xc,yc = proj(NX.array(cntrylons),NX.array(cntrylats))
- xc = xc.tolist(); yc = yc.tolist()
- xc2,yc2 = proj(NX.array(cntrylons2),NX.array(cntrylats))
- xc3,yc3 = proj(NX.array(cntrylons3),NX.array(cntrylats))
- xc2 = xc2.tolist(); yc2 = yc2.tolist()
- xc3 = xc3.tolist(); yc3 = yc3.tolist()
- segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(cntrysegind[:-1],cntrysegind[1:])]
- segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(cntrysegind[:-1],cntrysegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
- segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(cntrysegind[:-1],cntrysegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
- self.cntrysegs = segments+segments2+segments3
-
- # same as above for state segments.
- xc,yc = proj(NX.array(statelons),NX.array(statelats))
- xc = xc.tolist(); yc = yc.tolist()
- xc2,yc2 = proj(NX.array(statelons2),NX.array(statelats))
- xc3,yc3 = proj(NX.array(statelons3),NX.array(statelats))
- xc2 = xc2.tolist(); yc2 = yc2.tolist()
- xc3 = xc3.tolist(); yc3 = yc3.tolist()
- segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(statesegind[:-1],statesegind[1:])]
- segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(statesegind[:-1],statesegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
- segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(statesegind[:-1],statesegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
- self.statesegs = segments+segments2+segments3
-
- # same as above for river segments.
- xc,yc = proj(NX.array(riverlons),NX.array(riverlats))
- xc = xc.tolist(); yc = yc.tolist()
- xc2,yc2 = proj(NX.array(riverlons2),NX.array(riverlats))
- xc3,yc3 = proj(NX.array(riverlons3),NX.array(riverlats))
- xc2 = xc2.tolist(); yc2 = yc2.tolist()
- xc3 = xc3.tolist(); yc3 = yc3.tolist()
- segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(riversegind[:-1],riversegind[1:])]
- segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(riversegind[:-1],riversegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
- segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(riversegind[:-1],riversegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
- self.riversegs = segments+segments2+segments3
-
- # store coast polygons for filling.
+ self.area_thresh = area_thresh
+ # define map boundary polygon (in lat/lon coordinates)
+ self._boundarypoly = self._getmapboundary()
+ # read in coastline polygons, only keeping those that
+ # intersect map boundary polygon.
+ self.coastsegs, self.coastpolygontypes = self._readboundarydata('gshhs')
+ # same for countries, states, rivers.
+ self.cntrysegs, types = self._readboundarydata('countries')
+ self.statesegs, types = self._readboundarydata('states')
+ self.riversegs, types = self._readboundarydata('rivers')
+ # for coastlines, reformat as polygons.
self.coastpolygons = []
- coastpolygonsll = []
- self.coastpolygontypes = []
- if projection in ['merc','mill']:
- xsp,ysp = proj(0.,-89.9) # s. pole coordinates.
- xa,ya = proj(0.,-68.) # edge of antarctica.
- x0,y0 = proj(0.,0.)
- xm360,ym360 = proj(-360.,0.)
- x360,y360 = proj(360.,0.)
- x720,y720 = proj(720.,0.)
- for seg,segtype,segll in zip(self.coastsegs,self.coastsegtypes,self.coastsegsll):
- x = [lon for lon,lat in seg]
- y = [lat for lon,lat in seg]
- lons = [lon for lon,lat in segll]
- lats = [lat for lon,lat in segll]
- # the antarctic polygon is a nuisance, since it
- # spans all longitudes, it's not closed and will not be filled
- # without some projection dependant tweaking.
- if projection == 'cyl':
- if x[-1] == 0.000 and y[-1] < -68.: # close antarctica
- x.append(0.)
- y.append(-90.0000)
- x.insert(0,360.)
- y.insert(0,-90)
- lons.append(0.)
- lats.append(-90.)
- lons.insert(0,360.)
- lats.insert(0,-90.)
- if x[-1] == 360.000 and y[-1] < -68.:
- x.append(360.)
- y.append(-90)
- x.insert(0,720.)
- y.insert(0,-90)
- lons.append(360.)
- lats.append(-90.)
- lons.insert(0,720.)
- lats.insert(0,-90.)
- if x[-1] == -360.000 and y[-1] < -68.:
- x.append(-360.)
- y.append(-90)
- x.insert(0,0.)
- y.insert(0,-90)
- lons.append(-360.)
- lats.append(-90.)
- lons.insert(0,0.)
- lats.insert(0,-90.)
- elif projection in ['merc','mill']:
- if math.fabs(x[-1]-x0) < 1. and y[-1] < ya: # close antarctica
- x.append(x0)
- y.append(ysp)
- x.insert(0,x360)
- y.insert(0,ysp)
- lons.append(0.)
- lats.append(-90.)
- lons.insert(0,360.)
- lats.insert(0,-90.)
- if math.fabs(x[-1]-x360) < 1. and y[-1] < ya:
- x.append(x360)
- y.append(ysp)
- x.insert(0,x720)
- y.insert(0,ysp)
- lons.append(360.)
- lats.append(-90.)
- lons.insert(0,720.)
- lats.insert(0,-90.)
- if math.fabs(x[-1]-xm360) < 1. and y[-1] < ya:
- x.append(xm360)
- y.append(ysp)
- x.insert(0,x0)
- y.insert(0,ysp)
- lons.append(-360.)
- lats.append(-90.)
- lons.insert(0,0.)
- lats.insert(0,-90.)
+ for xy in self.coastsegs:
+ x = [x1 for x1,x2 in xy]
+ y = [x2 for x1,x2 in xy]
self.coastpolygons.append((x,y))
- coastpolygonsll.append((lons,lats))
- self.coastpolygontypes.append(segtype)
- # remove those segments/polygons that don't intersect map region.
- coastsegs = []
- coastsegtypes = []
- for seg,segtype in zip(self.coastsegs,self.coastsegtypes):
- if self._insidemap_seg(seg):
- coastsegs.append(seg)
- coastsegtypes.append(segtype)
- self.coastsegs = coastsegs
- self.coastsegtypes = coastsegtypes
- polygons = []
- polygonsll = []
- polygontypes = []
- for poly,polytype,polyll in zip(self.coastpolygons,self.coastpolygontypes,coastpolygonsll):
- if self._insidemap_poly(poly,polyll):
- polygons.append(poly)
- polygontypes.append(polytype)
- polygonsll.append(polyll)
- self.coastpolygons = polygons
- coastpolygonsll = polygonsll
- self.coastpolygontypes = polygontypes
- states = []; rivers = []; countries = []
- for seg in self.cntrysegs:
- if self._insidemap_seg(seg):
- countries.append(seg)
- for seg in self.statesegs:
- if self._insidemap_seg(seg):
- states.append(seg)
- for seg in self.riversegs:
- if self._insidemap_seg(seg):
- rivers.append(seg)
- self.statesegs = states
- self.riversegs = rivers
- self.cntryegs = countries
-
- # split up segments that go outside projection limb
- coastsegs = []
- coastsegtypes = []
- for seg,segtype in zip(self.coastsegs,self.coastsegtypes):
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
- i1,i2 = self._splitseg(xx,yy)
- if i1 and i2:
- for i,j in zip(i1,i2):
- segment = zip(xx[i:j].tolist(),yy[i:j].tolist())
- coastsegs.append(segment)
- coastsegtypes.append(segtype)
- else:
- coastsegs.append(seg)
- coastsegs.append(segtype)
- self.coastsegs = coastsegs
- self.coastsegtypes = coastsegtypes
- states = []
- for seg in self.statesegs:
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
- i1,i2 = self._splitseg(xx,yy)
- if i1 and i2:
- for i,j in zip(i1,i2):
- segment = zip(xx[i:j].tolist(),yy[i:j].tolist())
- states.append(segment)
- else:
- states.append(seg)
- self.statesegs = states
- countries = []
- for seg in self.cntrysegs:
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
- i1,i2 = self._splitseg(xx,yy)
- if i1 and i2:
- for i,j in zip(i1,i2):
- segment = zip(xx[i:j].tolist(),yy[i:j].tolist())
- countries.append(segment)
- else:
- countries.append(seg)
- self.cntrysegs = countries
- rivers = []
- for seg in self.riversegs:
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
- i1,i2 = self._splitseg(xx,yy)
- if i1 and i2:
- for i,j in zip(i1,i2):
- segment = zip(xx[i:j].tolist(),yy[i:j].tolist())
- rivers.append(segment)
- else:
- rivers.append(seg)
- self.riversegs = rivers
-
- # split coastline segments that jump across entire plot.
- coastsegs = []
- coastsegtypes = []
- for seg,segtype in zip(self.coastsegs,self.coastsegtypes):
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
- xd = (xx[1:]-xx[0:-1])**2
- yd = (yy[1:]-yy[0:-1])**2
- dist = NX.sqrt(xd+yd)
- split = dist > 5000000.
- if NX.sum(split) and self.projection not in ['merc','cyl','mill']:
- ind = (NX.compress(split,squeeze(split*NX.indices(xd.shape)))+1).tolist()
- iprev = 0
- ind.append(len(xd))
- for i in ind:
- coastsegs.append(zip(xx[iprev:i],yy[iprev:i]))
- coastsegtypes.append(segtype)
- iprev = i
- else:
- coastsegs.append(seg)
- coastsegtypes.append(segtype)
- self.coastsegs = coastsegs
- self.coastsegtypes = coastsegtypes
-
- # special treatment of coastline polygons for
- # geostationary, orthographic, sinusoidal, mollweide and robinson.
- # (polygon clipping along projection limb)
- if self.projection in ['ortho','geos']:
- if self.projection == 'ortho':
- lat_0 = math.radians(self.projparams['lat_0'])
- else:
- lat_0 = 0.
- lon_0d = self.projparams['lon_0']
- lon_0 = math.radians(lon_0d)
- if self.projection == 'ortho':
- rad = self.rmajor
- else:
- # quadratic mean radius of ellipsoid.
- rad = math.sqrt((3.*self._width**2 + self._height**2)/4.)
- del_s = 50.
- gc = pyproj.Geod(a=self.rmajor,b=self.rminor)
- coastpolygons = []
- coastpolygontypes = []
- for poly,polytype,polyll in zip(self.coastpolygons,self.coastpolygontypes,coastpolygonsll):
- x = poly[0]
- y = poly[1]
- lons = polyll[0]
- lats = polyll[1]
- mask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))
- # replace values in polygons that are over the horizon.
- xsave = False
- ysave = False
- if NX.sum(mask):
- i1,i2 = self._splitseg(x,y,mask=mask)
- # loop over segments of polygon that are outside projection limb.
- for i,j in zip(i1,i2):
- # if it's not the rest of the polygon ...
- if i and j != len(x):
- # compute distance and azimuth between projection center
- # and last point inside project limb.
- az1,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lons[i]),math.radians(lats[i]),radians=True)
- # also compute lat, lon of that great circle, plus back
- # azimuth.
- lon1,lat1,az=gc.fwd(lon_0,lat_0,az1,0.5*math.pi*rad,radians=True)
- # compute distance and azimuth between projection center
- # and next point inside projection limb.
- az2,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lons[j]),math.radians(lats[j]),radians=True)
- # also compute lat, lon of that great circle, plus back
- # azimuth.
- lon2,lat2,az=gc.fwd(lon_0,lat_0,az2,0.5*math.pi*rad,radians=True)
- # compute distance between those two points.
- az12,az21,dist = gc.inv(lon1,lon2,lat1,lat2,radians=True)
- # compute set of equally space points del_s meters apart
- # along great circle between those two points (the last
- # inside the projection limb and the next point inside the
- # the projection limb).
- npoints = int((dist+0.5*1000.*del_s)/(1000.*del_s))
- if npoints < 2: npoints=2
- lonlats = gc.npts(math.degrees(lon2),math.degrees(lat2),math.degrees(lon1),math.degrees(lat1),npoints)
- lonstmp=[math.degrees(lon2)];latstmp=[math.degrees(lat2)]
- for lon,lat in lonlats:
- lonstmp.append(lon); latstmp.append(lat)
- lonstmp.append(math.degrees(lon1)); latstmp.append(math.degrees(lat1))
- # convert that set of points to projection coordinates.
- # replace the points in the polygon which were outside
- # the projection limb.
- xx, yy = self(lonstmp, latstmp)
- xnew = x[i:j] + xx
- ynew = y[i:j] + yy
- coastpolygons.append((xnew,ynew))
- coastpolygontypes.append(polytype)
- elif i == 0:
- xsave = x[0:j]
- ysave = y[0:j]
- lats_save = lats[0:j]
- lons_save = lons[0:j]
- # it's the entire rest of the polygon ...
- elif j == len(x):
- xnew = x[i:j] + xsave
- ynew = y[i:j] + ysave
- lonsnew = lons[i:j] + lons_save
- latsnew = lats[i:j] + lats_save
- az1,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lonsnew[0]),math.radians(latsnew[0]),radians=True)
- lon1,lat1,az=gc.fwd(lon_0,lat_0,az1,0.5*math.pi*rad,radians=True)
- az2,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lonsnew[-1]),math.radians(latsnew[-1]),radians=True)
- lon2,lat2,az=gc.fwd(lon_0,lat_0,az2,0.5*math.pi*rad,radians=True)
- az12,az21,dist = gc.inv(lon2,lat2,lon1,lat1,radians=True)
- npoints = int((dist+0.5*1000.*del_s)/(1000.*del_s))
- if npoints < 2: npoints=2
- lonlats = gc.npts(math.degrees(lon2),math.degrees(lat2),math.degrees(lon1),math.degrees(lat1),npoints)
- lonstmp=[math.degrees(lon2)];latstmp=[math.degrees(lat2)]
- for lon,lat in lonlats:
- lonstmp.append(lon); latstmp.append(lat)
- lonstmp.append(math.degrees(lon1));latstmp.append(math.degrees(lat1))
- xx, yy = self(lonstmp, latstmp)
- xnew = xnew + xx
- ynew = ynew + yy
- coastpolygons.append((xnew,ynew))
- coastpolygontypes.append(polytype)
- else: # no part of polygon outside projection limb.
- coastpolygons.append(poly)
- coastpolygontypes.append(polytype)
- self.coastpolygons = coastpolygons
- self.coastpolygontypes = coastpolygontypes
- elif self.projection in ['moll','robin','sinu']:
- lon_0 = self.projparams['lon_0']
- coastpolygons=[]
- for poly,polytype,polyll in zip(self.coastpolygons,self.coastpolygontypes,coastpolygonsll):
- x = poly[0]
- y = poly[1]
- lons = polyll[0]
- lats = polyll[1]
- xn=[]
- yn=[]
- # antarctic segment goes from 360 back to 0
- # reorder to go from lon_0-180 to lon_0+180.
- if lats[-1] < -68.0:
- lons.reverse()
- lats.reverse()
- xx,yy = self(lons,lats)
- xx = NX.array(xx); yy = NX.array(yy)
- xdist = NX.fabs(xx[1:]-xx[0:-1])
- if max(xdist) > 1000000:
- nmin = NX.argmax(xdist)+1
- xnew = NX.zeros(len(xx),NX.Float64)
- ynew = NX.zeros(len(xx),NX.Float64)
- lonsnew = len(xx)*[0]
- latsnew = len(xx)*[0]
- xnew[0:len(xx)-nmin] = xx[nmin:]
- ynew[0:len(xx)-nmin] = yy[nmin:]
- xnew[len(xx)-nmin:] = xx[0:nmin]
- ynew[len(xx)-nmin:] = yy[0:nmin]
- lonsnew[0:len(xx)-nmin] = lons[nmin:]
- latsnew[0:len(xx)-nmin] = lats[nmin:]
- lonsnew[len(xx)-nmin:] = lons[0:nmin]
- latsnew[len(xx)-nmin:] = lats[0:nmin]
- x = xnew.tolist(); y = ynew.tolist()
- lons = lonsnew; lats = latsnew
- else:
- x.reverse()
- y.reverse()
- # close polygon (add lines along left and right edges down to S pole)
- for phi in NX.arange(-89.999,lats[0],0.1):
- xx,yy = self(lon_0-179.99,phi)
- xn.append(xx); yn.append(yy)
- xn = xn+x
- yn = yn+y
- for phi in NX.arange(lats[-1],-89.999,-0.1):
- xx,yy = self(lon_0+179.99,phi)
- xn.append(xx); yn.append(yy)
- # move points outside map to edge of map
- # along constant latitude.
- else:
- for x,y,lon,lat in zip(x,y,lons,lats):
- if lon > lon_0+180 or lon < lon_0-180:
- if lon >= lon_0+180: lon=lon_0+180.
- if lon <= lon_0-180: lon=lon_0-180.
- xx,yy = self(lon,lat)
- xn.append(xx); yn.append(yy)
- else:
- xn.append(x); yn.append(y)
- coastpolygons.append((xn,yn))
- self.coastpolygons = coastpolygons
-
def _splitseg(self,xx,yy,mask=None):
"""split segment up around missing values (outside projection limb)"""
if mask is None:
@@ -1266,6 +770,138 @@
"""
return self.projtran.makegrid(nx,ny,returnxy=returnxy)
+ def _readboundarydata(self,name):
+ msg = """
+Unable to open boundary dataset file. Only the 'crude', 'low'
+and 'intermediate' resolution datasets are installed by default. If you
+are requesting a 'high' resolution dataset, you need to download
+and install those files manually (see the basemap README for details)."""
+ try:
+ bdatfile = open(os.path.join(basemap_datadir,name+'_'+self.resolution+'.dat'),'rb')
+ bdatmetafile = open(os.path.join(basemap_datadir,name+'meta_'+self.resolution+'.dat'),'r')
+ except:
+ raise IOError, msg
+ polygons = []
+ polygon_types = []
+ for line in bdatmetafile:
+ linesplit = line.split()
+ area = float(linesplit[1])
+ type = int(linesplit[0])
+ south = float(linesplit[2])
+ north = float(linesplit[3])
+ if area < 0.: area = 1.e30
+ useit = self.latmax>=south and self.latmin<=north and area>self.area_thresh
+ # skip Antartica for now.
+ if name == 'gshhs' and south < -60.: useit=False
+ if useit:
+ offsetbytes = int(linesplit[4])
+ bytecount = int(linesplit[5])
+ bdatfile.seek(offsetbytes,0)
+ polystring = bdatfile.read(bytecount)
+ poly = wkb.loads(polystring)
+ if poly.intersects(self._boundarypoly):
+ #a = npy.asarray(self._boundarypoly.boundary)
+ #b = npy.asarray(poly.boundary)
+ #import pylab
+ #pylab.fill(a[:,0],a[:,1],'r')
+ #pylab.fill(b[:,0],b[:,1],'b')
+ #pylab.show()
+ poly = poly.intersection(self._boundarypoly)
+ if hasattr(poly,'geoms'):
+ geoms = poly.geoms
+ else:
+ geoms = [poly]
+ for psub in geoms:
+ if name == 'gshhs':
+ b = npy.asarray(psub.boundary)
+ else:
+ b = npy.asarray(psub.coords)
+ blons = b[:,0]; blats = b[:,1]
+ bx, by = self(blons, blats)
+ #if (bx > 1.20).any() or (by > 1.e20).any():
+ # continue
+ polygons.append(zip(bx,by))
+ polygon_types.append(type)
+ return polygons, polygon_types
+
+
+ def _getmapboundary(self):
+ """
+ define map boundary polygon (in lat/lon coordinates)
+ """
+ dtheta = 0.1
+ dx = (self.xmax-self.xmin)/100.
+ dy = (self.ymax-self.ymin)/100.
+ if self.projection == 'ortho' and self._fulldisk:
+ # circular region.
+ thetas = npy.arange(0.,2.*npy.pi,dtheta)
+ radius = self.rmajor
+ x = radius*npy.cos(thetas) + 0.5*self.xmax
+ y = radius*npy.sin(thetas) + 0.5*self.ymax
+ elif self.projection == 'geos' and self._fulldisk:
+ # elliptical region
+ thetas = npy.arange(0.,2.*npy.pi+0.5*dtheta,dtheta)
+ rminor = self._height
+ rmajor = self._width
+ x = rmajor*npy.cos(thetas) + 0.5*self.xmax
+ y = rminor*npy.sin(thetas) + 0.5*self.ymax
+ elif self.projection in ['moll','robin','sinu']:
+ # quasi-elliptical region.
+ x = []; y = []
+ # left side
+ lats = NX.arange(-89.9,89.9+dtheta,dtheta).tolist()
+ lons = len(lats)*[self.projparams['lon_0']-179.9]
+ x,y = self(lons,lats)
+ # top.
+ lons = NX.arange(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179+dtheta,dtheta).tolist()
+ lats = len(lons)*[89.9]
+ xx,yy = self(lons,lats)
+ x = x+xx; y = y+yy
+ # right side
+ lats = NX.arange(89.9,-89.9-dtheta,-dtheta).tolist()
+ lons = len(lats)*[self.projparams['lon_0']+179.9]
+ xx,yy = self(lons,lats)
+ x = x+xx; y = y+yy
+ # bottom.
+ lons = NX.arange(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180-dtheta,-dtheta).tolist()
+ lats = len(lons)*[-89.9]
+ xx,yy = self(lons,lats)
+ x = x+xx; y = y+yy
+ x = npy.array(x,npy.float64)
+ y = npy.array(y,npy.float64)
+ else: # all other projections are rectangular.
+ # left side (x = xmin, ymin <= y <= ymax)
+ yy = npy.arange(self.ymin, self.ymax+0.5*dy, dy)
+ x = len(yy)*[self.xmin]; y = yy.tolist()
+ # top (y = ymax, xmin <= x <= xmax)
+ xx = npy.arange(self.xmin, self.xmax+0.5*dx, dx)
+ x = x + xx.tolist()
+ y = y + len(xx)*[self.ymax]
+ # right side (x = xmax, ymin <= y <= ymax)
+ yy = npy.arange(self.ymax, self.ymin-0.5*dy, -dy)
+ x = x + len(yy)*[self.xmax]; y = y + yy.tolist()
+ # bottom (y = ymin, xmin <= x <= xmax)
+ xx = npy.arange(self.xmax, self.xmin-0.5*dx, -dx)
+ x = x + xx.tolist()
+ y = y + len(xx)*[self.ymin]
+ x = npy.array(x,npy.float64)
+ y = npy.array(y,npy.float64)
+ lons, lats = self(x,y,inverse=True)
+ # fix lons so there are no jumps.
+ n = 1
+ lonprev = lons[0]
+ for lon,lat in zip(lons[1:],lats[1:]):
+ if npy.abs(lon-lonprev) > 90.:
+ if lonprev < 0:
+ lon = lon - 360.
+ else:
+ lon = lon + 360
+ lons[n] = lon
+ lonprev = lon
+ n = n + 1
+ return PolygonShape(zip(lons,lats))
+
+
def drawmapboundary(self,color='k',linewidth=1.0,ax=None):
"""
draw boundary around map projection region. If ax=None (default),
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-01 03:30:11
|
Revision: 4079
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4079&view=rev
Author: jswhit
Date: 2007-10-31 20:30:07 -0700 (Wed, 31 Oct 2007)
Log Message:
-----------
delete uneeded methods
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-01 03:26:12 UTC (rev 4078)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-01 03:30:07 UTC (rev 4079)
@@ -688,59 +688,14 @@
self.cntrysegs, types = self._readboundarydata('countries')
self.statesegs, types = self._readboundarydata('states')
self.riversegs, types = self._readboundarydata('rivers')
- # for coastlines, reformat as polygons.
+ # for coastlines, reformat for use in
+ # matplotlib.patches.Polygon.
self.coastpolygons = []
for xy in self.coastsegs:
x = [x1 for x1,x2 in xy]
y = [x2 for x1,x2 in xy]
self.coastpolygons.append((x,y))
- def _splitseg(self,xx,yy,mask=None):
- """split segment up around missing values (outside projection limb)"""
- if mask is None:
- mask = (NX.logical_or(NX.greater_equal(xx,1.e20),NX.greater_equal(yy,1.e20))).tolist()
- i1=[]; i2=[]
- mprev = 1
- for i,m in enumerate(mask):
- if not m and mprev:
- i1.append(i)
- if m and not mprev:
- i2.append(i)
- mprev = m
- if not mprev: i2.append(len(mask))
- if len(i1) != len(i2):
- raise ValueError,'error in splitting boundary segments'
- return i1,i2
-
- def _insidemap_seg(self,seg):
- """returns True if any point in segment is inside map region"""
- xx = [x for x,y in seg]
- yy = [y for x,y in seg]
- isin = False
- for x,y in zip(xx,yy):
- if x >= self.xmin and x <= self.xmax and y >= self.ymin and y <= self.ymax:
- isin = True
- break
- return isin
-
- def _insidemap_poly(self,poly,polyll):
- """returns True if any point in polygon is inside map region"""
- isin = False
- xx = poly[0]; yy = poly[1]
- if self.projection in ['moll','robin','sinu']:
- lon_0 = self.projparams['lon_0']
- lons = polyll[0]
- for lon in lons:
- if lon < lon_0+180 and lon > lon_0-180:
- isin = True
- break
- else:
- for x,y in zip(xx,yy):
- if x >= self.xmin and x <= self.xmax and y >= self.ymin and y <= self.ymax:
- isin = True
- break
- return isin
-
def __call__(self,x,y,inverse=False):
"""
Calling a Basemap class instance with the arguments lon, lat will
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-01 12:05:31
|
Revision: 4082
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4082&view=rev
Author: jswhit
Date: 2007-11-01 05:05:30 -0700 (Thu, 01 Nov 2007)
Log Message:
-----------
fix antarctica
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-01 07:21:00 UTC (rev 4081)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-01 12:05:30 UTC (rev 4082)
@@ -747,20 +747,36 @@
if area < 0.: area = 1.e30
useit = self.latmax>=south and self.latmin<=north and area>self.area_thresh
# skip Antartica for now.
- if name == 'gshhs' and south < -60.: useit=False
if useit:
offsetbytes = int(linesplit[4])
bytecount = int(linesplit[5])
bdatfile.seek(offsetbytes,0)
polystring = bdatfile.read(bytecount)
poly = wkb.loads(polystring)
+ # close Antarctica for cylindrical projections
+
+ if name == 'gshhs' and self.projection in ['cyl','merc','mill']:
+ b = npy.asarray(poly.boundary)
+ lons = b[:,0].tolist()
+ lats = b[:,1].tolist()
+ if (math.fabs(lons[0]-360.) < 1.e-5 or
+ math.fabs(lons[0]+360.) or
+ math.fabs(lons[0]-0.) < 1.e-5) and lats[-1] < -68.:
+ lons = lons[:-2]
+ lats = lats[:-2]
+ lonstart,latstart = lons[0], lats[0]
+ lonend,latend = lons[-1], lats[-1]
+ lons.insert(0,lonstart)
+ lats.insert(0,-90.)
+ lons.append(lonend)
+ lats.append(-90.)
+ if south < -68:
+ poly = PolygonShape(zip(lons,lats))
+ #b = npy.asarray(poly.boundary)
+ #import pylab
+ #pylab.fill(b[:,0],b[:,1],'b')
+ #pylab.show()
if poly.intersects(self._boundarypoly):
- #a = npy.asarray(self._boundarypoly.boundary)
- #b = npy.asarray(poly.boundary)
- #import pylab
- #pylab.fill(a[:,0],a[:,1],'r')
- #pylab.fill(b[:,0],b[:,1],'b')
- #pylab.show()
poly = poly.intersection(self._boundarypoly)
if hasattr(poly,'geoms'):
geoms = poly.geoms
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-01 12:10:45
|
Revision: 4083
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4083&view=rev
Author: jswhit
Date: 2007-11-01 05:10:40 -0700 (Thu, 01 Nov 2007)
Log Message:
-----------
partial fix for Antarctica.
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-01 12:05:30 UTC (rev 4082)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-01 12:10:40 UTC (rev 4083)
@@ -770,7 +770,6 @@
lats.insert(0,-90.)
lons.append(lonend)
lats.append(-90.)
- if south < -68:
poly = PolygonShape(zip(lons,lats))
#b = npy.asarray(poly.boundary)
#import pylab
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-01 14:02:34
|
Revision: 4086
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4086&view=rev
Author: jswhit
Date: 2007-11-01 07:02:29 -0700 (Thu, 01 Nov 2007)
Log Message:
-----------
more fixes for Antarctica
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-01 13:01:36 UTC (rev 4085)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-01 14:02:29 UTC (rev 4086)
@@ -757,24 +757,29 @@
if name == 'gshhs' and self.projection in ['cyl','merc','mill']:
b = npy.asarray(poly.boundary)
- lons = b[:,0].tolist()
- lats = b[:,1].tolist()
- if (math.fabs(lons[0]-360.) < 1.e-5 or
- math.fabs(lons[0]+360.) or
- math.fabs(lons[0]-0.) < 1.e-5) and lats[-1] < -68.:
- lons = lons[:-2]
- lats = lats[:-2]
- lonstart,latstart = lons[0], lats[0]
- lonend,latend = lons[-1], lats[-1]
- lons.insert(0,lonstart)
- lats.insert(0,-90.)
- lons.append(lonend)
- lats.append(-90.)
- poly = PolygonShape(zip(lons,lats))
- #b = npy.asarray(poly.boundary)
- #import pylab
- #pylab.fill(b[:,0],b[:,1],'b')
- #pylab.show()
+ lons = b[:,0]
+ lats = b[:,1]
+ if south < -68:
+ if math.fabs(lons[0]+0.) < 1.e-5:
+ lons1 = lons[:-2][::-1]
+ lats1 = lats[:-2][::-1]
+ lons2 = lons1 + 360.
+ lons3 = lons2 + 360.
+ lons = lons1.tolist()+lons2.tolist()+lons3.tolist()
+ lats = lats1.tolist()+lats1.tolist()+lats1.tolist()
+ lonstart,latstart = lons[0], lats[0]
+ lonend,latend = lons[-1], lats[-1]
+ lons.insert(0,lonstart)
+ lats.insert(0,-90.)
+ lons.append(lonend)
+ lats.append(-90.)
+ poly = PolygonShape(zip(lons,lats))
+ #b = npy.asarray(poly.boundary)
+ #import pylab
+ #pylab.fill(b[:,0],b[:,1],'b')
+ #pylab.show()
+ else:
+ continue
if poly.intersects(self._boundarypoly):
poly = poly.intersection(self._boundarypoly)
if hasattr(poly,'geoms'):
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-01 16:25:00
|
Revision: 4091
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4091&view=rev
Author: jswhit
Date: 2007-11-01 09:24:53 -0700 (Thu, 01 Nov 2007)
Log Message:
-----------
cleanup
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-01 15:17:24 UTC (rev 4090)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-01 16:24:53 UTC (rev 4091)
@@ -747,14 +747,14 @@
if area < 0.: area = 1.e30
useit = self.latmax>=south and self.latmin<=north and area>self.area_thresh
# skip Antartica for now.
+ #if south < -68: useit = False
if useit:
offsetbytes = int(linesplit[4])
bytecount = int(linesplit[5])
bdatfile.seek(offsetbytes,0)
polystring = bdatfile.read(bytecount)
poly = wkb.loads(polystring)
- # close Antarctica for cylindrical projections
-
+ # close Antartica for cylindrical projections.
if name == 'gshhs' and self.projection in ['cyl','merc','mill']:
b = npy.asarray(poly.boundary)
lons = b[:,0]
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-02 23:21:25
|
Revision: 4098
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4098&view=rev
Author: jswhit
Date: 2007-11-02 16:21:06 -0700 (Fri, 02 Nov 2007)
Log Message:
-----------
fixed for map projection regions that are not polygons in lat/lon coords
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-02 18:45:38 UTC (rev 4097)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-02 23:21:06 UTC (rev 4098)
@@ -10,7 +10,6 @@
from matplotlib.lines import Line2D
import pyproj, sys, os, math, dbflib
from proj import Proj
-import matplotlib.numerix as NX
from matplotlib.numerix import ma
import numpy as npy
from numpy import linspace
@@ -18,17 +17,15 @@
from matplotlib.cbook import popd, is_scalar
from shapelib import ShapeFile
from shapely.geometry import Polygon as PolygonShape
-from shapely import wkb
+from shapely.geometry import LineString as LineShape
+from shapely.geometry import Point as PointShape
+from shapely import wkb, wkt
# basemap data files now installed in lib/matplotlib/toolkits/basemap/data
basemap_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
__version__ = '0.9.7'
-# test to see numerix set to use numpy (if not, raise an error)
-if NX.which[0] != 'numpy':
- raise ImportError("numerix must be set to numpy")
-
class Basemap(object):
"""
@@ -655,8 +652,8 @@
self.latmax = 90.
else:
lons, lats = self.makegrid(101,101)
- self.latmin = min(NX.ravel(lats))
- self.latmax = max(NX.ravel(lats))
+ self.latmin = lats.min()
+ self.latmax = lats.max()
# if ax == None, pylab.gca may be used.
self.ax = ax
@@ -680,7 +677,7 @@
raise ValueError, "boundary resolution must be one of 'c','l','i' or 'h'"
self.area_thresh = area_thresh
# define map boundary polygon (in lat/lon coordinates)
- self._boundarypoly = self._getmapboundary()
+ self._boundarypolyll, self._boundarypolyxy = self._getmapboundary()
# read in coastline polygons, only keeping those that
# intersect map boundary polygon.
self.coastsegs, self.coastpolygontypes = self._readboundarydata('gshhs')
@@ -738,6 +735,13 @@
raise IOError, msg
polygons = []
polygon_types = []
+ # see if map projection region polygon contains a pole.
+ NPole = PointShape(self(0.,90.))
+ SPole = PointShape(self(0.,-90.))
+ hasNP = self._boundarypolyxy.contains(NPole)
+ hasSP = self._boundarypolyxy.contains(SPole)
+ containsPole = hasNP or hasSP
+ # iterate over boundary geometries.
for line in bdatmetafile:
linesplit = line.split()
area = float(linesplit[1])
@@ -746,57 +750,117 @@
north = float(linesplit[3])
if area < 0.: area = 1.e30
useit = self.latmax>=south and self.latmin<=north and area>self.area_thresh
- # skip Antartica for now.
- #if south < -68: useit = False
if useit:
offsetbytes = int(linesplit[4])
bytecount = int(linesplit[5])
bdatfile.seek(offsetbytes,0)
polystring = bdatfile.read(bytecount)
poly = wkb.loads(polystring)
- # close Antartica for cylindrical projections.
- if name == 'gshhs' and self.projection in ['cyl','merc','mill']:
- b = npy.asarray(poly.boundary)
- lons = b[:,0]
- lats = b[:,1]
- if south < -68:
- if math.fabs(lons[0]+0.) < 1.e-5:
- lons1 = lons[:-2][::-1]
- lats1 = lats[:-2][::-1]
- lons2 = lons1 + 360.
- lons3 = lons2 + 360.
- lons = lons1.tolist()+lons2.tolist()+lons3.tolist()
- lats = lats1.tolist()+lats1.tolist()+lats1.tolist()
- lonstart,latstart = lons[0], lats[0]
- lonend,latend = lons[-1], lats[-1]
- lons.insert(0,lonstart)
- lats.insert(0,-90.)
- lons.append(lonend)
- lats.append(-90.)
- poly = PolygonShape(zip(lons,lats))
- #b = npy.asarray(poly.boundary)
- #import pylab
- #pylab.fill(b[:,0],b[:,1],'b')
- #pylab.show()
+ # if map boundary polygon is a valid one in lat/lon
+ # coordinates (i.e. it does not contain either pole),
+ # the intersections of the boundary geometries
+ # and the map projection region can be computed before
+ # transforming the boundary geometry to map projection
+ # coordinates (this saves time, especially for small map
+ # regions and high-resolution boundary geometries).
+ if not containsPole:
+ # close Antarctica for cylindrical projections.
+ if name == 'gshhs' and self.projection in ['cyl','merc','mill']:
+ b = npy.asarray(poly.boundary)
+ lons = b[:,0]
+ lats = b[:,1]
+ if south < -68:
+ if math.fabs(lons[0]+0.) < 1.e-5:
+ lons1 = lons[:-2][::-1]
+ lats1 = lats[:-2][::-1]
+ lons2 = lons1 + 360.
+ lons3 = lons2 + 360.
+ lons = lons1.tolist()+lons2.tolist()+lons3.tolist()
+ lats = lats1.tolist()+lats1.tolist()+lats1.tolist()
+ lonstart,latstart = lons[0], lats[0]
+ lonend,latend = lons[-1], lats[-1]
+ lons.insert(0,lonstart)
+ lats.insert(0,-90.)
+ lons.append(lonend)
+ lats.append(-90.)
+ poly = PolygonShape(zip(lons,lats))
+ #b = npy.asarray(poly.boundary)
+ #import pylab
+ #pylab.fill(b[:,0],b[:,1],'b')
+ #pylab.show()
+ else:
+ continue
+ if poly.intersects(self._boundarypolyll):
+ poly = poly.intersection(self._boundarypolyll)
+ if hasattr(poly,'geoms'):
+ geoms = poly.geoms
else:
- continue
- if poly.intersects(self._boundarypoly):
- poly = poly.intersection(self._boundarypoly)
- if hasattr(poly,'geoms'):
- geoms = poly.geoms
+ geoms = [poly]
+ for psub in geoms:
+ if name == 'gshhs':
+ b = npy.asarray(psub.boundary)
+ else:
+ b = npy.asarray(psub.coords)
+ blons = b[:,0]; blats = b[:,1]
+ bx, by = self(blons, blats)
+ #if (bx > 1.20).any() or (by > 1.e20).any():
+ # continue
+ polygons.append(zip(bx,by))
+ polygon_types.append(type)
+ # if map boundary polygon is not valid in lat/lon
+ # coordinates, compute intersection between map
+ # projection region and boundary geometries in map
+ # projection coordinates.
+ else:
+ if name == 'gshhs':
+ b = npy.asarray(poly.boundary)
else:
- geoms = [poly]
- for psub in geoms:
- if name == 'gshhs':
- b = npy.asarray(psub.boundary)
+ b = npy.asarray(poly.coords)
+ bx, by = self(b[:,0], b[:,1])
+ mask = npy.logical_and(bx<1.e20,by<1.e20)
+ if sum(mask) <= 1: continue
+ if name == 'gshhs':
+ poly = PolygonShape(zip(bx,by))
+ else:
+ # remove parts of geometry that are undefined
+ # in this map projection.
+ bx = npy.compress(mask, bx)
+ by = npy.compress(mask, by)
+ poly = LineShape(zip(bx,by))
+ if 0:
+ import pylab
+ print poly
+ pp = self._boundarypolyxy.intersection(poly)
+ print name, pp
+ a = npy.asarray(self._boundarypolyxy.boundary)
+ pylab.plot(a[:,0],a[:,1],'b')
+ pylab.plot(bx,by,'r')
+ if hasattr(pp,'geoms'):
+ px = pp.geoms
else:
- b = npy.asarray(psub.coords)
- blons = b[:,0]; blats = b[:,1]
- bx, by = self(blons, blats)
- #if (bx > 1.20).any() or (by > 1.e20).any():
- # continue
- polygons.append(zip(bx,by))
- polygon_types.append(type)
+ px = [pp]
+ for p in px:
+ c = npy.asarray(p.boundary)
+ pylab.fill(c[:,0],c[:,1],'g')
+ pylab.show()
+ raise SystemExit(0)
+ if self._boundarypolyxy.intersects(poly):
+ try:
+ poly = self._boundarypolyxy.intersection(poly)
+ except:
+ continue
+ if hasattr(poly,'geoms'):
+ geoms = poly.geoms
+ else:
+ geoms = [poly]
+ for psub in geoms:
+ if name == 'gshhs':
+ b = npy.asarray(psub.boundary)
+ else:
+ b = npy.asarray(psub.coords)
+ polygons.append(zip(b[:,0],b[:,1]))
+ polygon_types.append(type)
+ print name,len(polygons)
return polygons, polygon_types
@@ -804,77 +868,89 @@
"""
define map boundary polygon (in lat/lon coordinates)
"""
- dtheta = 0.1
- dx = (self.xmax-self.xmin)/100.
- dy = (self.ymax-self.ymin)/100.
+ nx = 100
+ ny = 100
if self.projection == 'ortho' and self._fulldisk:
# circular region.
- thetas = npy.arange(0.,2.*npy.pi,dtheta)
+ thetas = linspace(0.,2.*npy.pi,2*nx*ny)[:-1]
radius = self.rmajor
x = radius*npy.cos(thetas) + 0.5*self.xmax
y = radius*npy.sin(thetas) + 0.5*self.ymax
+ boundaryxy = PolygonShape(zip(x,y))
elif self.projection == 'geos' and self._fulldisk:
# elliptical region
- thetas = npy.arange(0.,2.*npy.pi+0.5*dtheta,dtheta)
+ thetas = linspace(0.,2.*npy.pi,2*nx*ny)[:-1]
rminor = self._height
rmajor = self._width
x = rmajor*npy.cos(thetas) + 0.5*self.xmax
y = rminor*npy.sin(thetas) + 0.5*self.ymax
+ boundaryxy = PolygonShape(zip(x,y))
elif self.projection in ['moll','robin','sinu']:
# quasi-elliptical region.
x = []; y = []
# left side
- lats = NX.arange(-89.9,89.9+dtheta,dtheta).tolist()
+ lats = linspace(-89.9,89.9,ny).tolist()
lons = len(lats)*[self.projparams['lon_0']-179.9]
x,y = self(lons,lats)
# top.
- lons = NX.arange(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179+dtheta,dtheta).tolist()
+ lons = linspace(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179,nx).tolist()
lats = len(lons)*[89.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
# right side
- lats = NX.arange(89.9,-89.9-dtheta,-dtheta).tolist()
+ lats = linspace(89.9,-89.9,ny).tolist()
lons = len(lats)*[self.projparams['lon_0']+179.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
# bottom.
- lons = NX.arange(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180-dtheta,-dtheta).tolist()
+ lons = linspace(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180).tolist()
lats = len(lons)*[-89.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
x = npy.array(x,npy.float64)
y = npy.array(y,npy.float64)
+ boundaryxy = PolygonShape(zip(x,y))
else: # all other projections are rectangular.
# left side (x = xmin, ymin <= y <= ymax)
- yy = npy.arange(self.ymin, self.ymax+0.5*dy, dy)
+ yy = linspace(self.ymin, self.ymax, ny)[:-1]
x = len(yy)*[self.xmin]; y = yy.tolist()
# top (y = ymax, xmin <= x <= xmax)
- xx = npy.arange(self.xmin, self.xmax+0.5*dx, dx)
+ xx = npy.linspace(self.xmin, self.xmax, nx)[:-1]
x = x + xx.tolist()
y = y + len(xx)*[self.ymax]
# right side (x = xmax, ymin <= y <= ymax)
- yy = npy.arange(self.ymax, self.ymin-0.5*dy, -dy)
+ yy = linspace(self.ymax, self.ymin, ny)[:-1]
x = x + len(yy)*[self.xmax]; y = y + yy.tolist()
# bottom (y = ymin, xmin <= x <= xmax)
- xx = npy.arange(self.xmax, self.xmin-0.5*dx, -dx)
+ xx = linspace(self.xmax, self.xmin, nx)[:-1]
x = x + xx.tolist()
y = y + len(xx)*[self.ymin]
x = npy.array(x,npy.float64)
y = npy.array(y,npy.float64)
- lons, lats = self(x,y,inverse=True)
- # fix lons so there are no jumps.
- n = 1
- lonprev = lons[0]
- for lon,lat in zip(lons[1:],lats[1:]):
- if npy.abs(lon-lonprev) > 90.:
- if lonprev < 0:
- lon = lon - 360.
- else:
- lon = lon + 360
- lons[n] = lon
- lonprev = lon
- n = n + 1
- return PolygonShape(zip(lons,lats))
+ boundaryxy = PolygonShape([(self.xmin,self.ymin),\
+ (self.xmin,self.ymax),\
+ (self.xmax,self.ymax),\
+ (self.xmax,self.ymin)])
+ if self.projection in ['mill','merc','cyl']:
+ lons = [self.llcrnrlon, self.llcrnrlon, self.urcrnrlon, self.urcrnrlon]
+ lats = [self.llcrnrlat, self.urcrnrlat, self.urcrnrlat, self.llcrnrlat]
+ else:
+ lons, lats = self(x,y,inverse=True)
+ # fix lons so there are no jumps.
+ n = 1
+ lonprev = lons[0]
+ for lon,lat in zip(lons[1:],lats[1:]):
+ if npy.abs(lon-lonprev) > 90.:
+ if lonprev < 0:
+ lon = lon - 360.
+ else:
+ lon = lon + 360
+ lons[n] = lon
+ lonprev = lon
+ n = n + 1
+ boundaryll = PolygonShape(zip(lons,lats))
+ #print 'map projection region',boundaryll.geom_type,boundaryll.is_valid
+ return PolygonShape(zip(lons,lats)), boundaryxy
def drawmapboundary(self,color='k',linewidth=1.0,ax=None):
@@ -912,21 +988,21 @@
ellps.set_clip_on(False)
elif self.projection in ['moll','robin','sinu']: # elliptical region.
# left side
- lats = NX.arange(-89.9,89.9+dtheta,dtheta).tolist()
+ lats = npy.arange(-89.9,89.9+dtheta,dtheta).tolist()
lons = len(lats)*[self.projparams['lon_0']-179.9]
x,y = self(lons,lats)
# top.
- lons = NX.arange(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179+dtheta,dtheta).tolist()
+ lons = npy.arange(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179+dtheta,dtheta).tolist()
lats = len(lons)*[89.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
# right side
- lats = NX.arange(89.9,-89.9-dtheta,-dtheta).tolist()
+ lats = npy.arange(89.9,-89.9-dtheta,-dtheta).tolist()
lons = len(lats)*[self.projparams['lon_0']+179.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
# bottom.
- lons = NX.arange(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180-dtheta,-dtheta).tolist()
+ lons = npy.arange(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180-dtheta,-dtheta).tolist()
lats = len(lons)*[-89.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
@@ -970,18 +1046,18 @@
axisbgc = ax.get_axis_bgcolor()
np = 0
for x,y in self.coastpolygons:
- xa = NX.array(x,NX.Float32)
- ya = NX.array(y,NX.Float32)
+ xa = npy.array(x,npy.float32)
+ ya = npy.array(y,npy.float32)
# check to see if all four corners of domain in polygon (if so,
# don't draw since it will just fill in the whole map).
delx = 10; dely = 10
if self.projection in ['cyl']:
delx = 0.1
dely = 0.1
- test1 = NX.fabs(xa-self.urcrnrx) < delx
- test2 = NX.fabs(xa-self.llcrnrx) < delx
- test3 = NX.fabs(ya-self.urcrnry) < dely
- test4 = NX.fabs(ya-self.llcrnry) < dely
+ test1 = npy.fabs(xa-self.urcrnrx) < delx
+ 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)
@@ -1277,9 +1353,9 @@
xoffset = (self.urcrnrx-self.llcrnrx)/100.
if self.projection in ['merc','cyl','mill','moll','robin','sinu']:
- lons = NX.arange(self.llcrnrlon,self.urcrnrlon+0.01,0.01)
+ lons = npy.arange(self.llcrnrlon,self.urcrnrlon+0.01,0.01)
else:
- lons = NX.arange(0,360.01,0.01)
+ lons = npy.arange(0,360.01,0.01)
# make sure latmax degree parallel is drawn if projection not merc or cyl or miller
try:
circlesl = circles.tolist()
@@ -1293,24 +1369,24 @@
xdelta = 0.01*(self.xmax-self.xmin)
ydelta = 0.01*(self.ymax-self.ymin)
for circ in circlesl:
- lats = circ*NX.ones(len(lons),NX.Float32)
+ lats = circ*npy.ones(len(lons),npy.float32)
x,y = self(lons,lats)
# remove points outside domain.
- testx = NX.logical_and(x>=self.xmin-xdelta,x<=self.xmax+xdelta)
- x = NX.compress(testx, x)
- y = NX.compress(testx, y)
- testy = NX.logical_and(y>=self.ymin-ydelta,y<=self.ymax+ydelta)
- x = NX.compress(testy, x)
- y = NX.compress(testy, y)
+ testx = npy.logical_and(x>=self.xmin-xdelta,x<=self.xmax+xdelta)
+ x = npy.compress(testx, x)
+ y = npy.compress(testx, y)
+ testy = npy.logical_and(y>=self.ymin-ydelta,y<=self.ymax+ydelta)
+ x = npy.compress(testy, x)
+ y = npy.compress(testy, y)
if len(x) > 1 and len(y) > 1:
# split into separate line segments if necessary.
# (not necessary for mercator or cylindrical or miller).
xd = (x[1:]-x[0:-1])**2
yd = (y[1:]-y[0:-1])**2
- dist = NX.sqrt(xd+yd)
+ dist = npy.sqrt(xd+yd)
split = dist > 500000.
- if NX.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
- ind = (NX.compress(split,squeeze(split*NX.indices(xd.shape)))+1).tolist()
+ if npy.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
+ ind = (npy.compress(split,squeeze(split*npy.indices(xd.shape)))+1).tolist()
xl = []
yl = []
iprev = 0
@@ -1355,10 +1431,10 @@
if self.projection == 'moll' and yy[0] < 1.e-4:
yy[0] = 1.e-4
if side == 'l':
- lons,lats = self(self.llcrnrx*NX.ones(yy.shape,NX.Float32),yy,inverse=True)
+ lons,lats = self(self.llcrnrx*npy.ones(yy.shape,npy.float32),yy,inverse=True)
lons = lons.tolist(); lats = lats.tolist()
else:
- lons,lats = self(self.urcrnrx*NX.ones(yy.shape,NX.Float32),yy,inverse=True)
+ lons,lats = self(self.urcrnrx*npy.ones(yy.shape,npy.float32),yy,inverse=True)
lons = lons.tolist(); lats = lats.tolist()
if max(lons) > 1.e20 or max(lats) > 1.e20:
raise ValueError,'inverse transformation undefined - please adjust the map projection region'
@@ -1368,10 +1444,10 @@
nmax = int((self.xmax-self.xmin)/dx+1)
xx = linspace(self.llcrnrx,self.urcrnrx,nmax)
if side == 'b':
- lons,lats = self(xx,self.llcrnry*NX.ones(xx.shape,NX.Float32),inverse=True)
+ lons,lats = self(xx,self.llcrnry*npy.ones(xx.shape,npy.float32),inverse=True)
lons = lons.tolist(); lats = lats.tolist()
else:
- lons,lats = self(xx,self.urcrnry*NX.ones(xx.shape,NX.Float32),inverse=True)
+ lons,lats = self(xx,self.urcrnry*npy.ones(xx.shape,npy.float32),inverse=True)
lons = lons.tolist(); lats = lats.tolist()
if max(lons) > 1.e20 or max(lats) > 1.e20:
raise ValueError,'inverse transformation undefined - please adjust the map projection region'
@@ -1394,7 +1470,7 @@
latlabstr = u'-%s\N{DEGREE SIGN}'%fmt
else:
latlabstr = u'%s\N{DEGREE SIGN}S'%fmt
- latlab = latlabstr%NX.fabs(lat)
+ latlab = latlabstr%npy.fabs(lat)
elif lat>0:
if rcParams['text.usetex']:
if labelstyle=='+/-':
@@ -1493,30 +1569,30 @@
xoffset = (self.urcrnrx-self.llcrnrx)/100.
if self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
- lats = NX.arange(-latmax,latmax+0.01,0.01)
+ lats = npy.arange(-latmax,latmax+0.01,0.01)
else:
- lats = NX.arange(-90,90.01,0.01)
+ lats = npy.arange(-90,90.01,0.01)
xdelta = 0.01*(self.xmax-self.xmin)
ydelta = 0.01*(self.ymax-self.ymin)
for merid in meridians:
- lons = merid*NX.ones(len(lats),NX.Float32)
+ lons = merid*npy.ones(len(lats),npy.float32)
x,y = self(lons,lats)
# remove points outside domain.
- testx = NX.logical_and(x>=self.xmin-xdelta,x<=self.xmax+xdelta)
- x = NX.compress(testx, x)
- y = NX.compress(testx, y)
- testy = NX.logical_and(y>=self.ymin-ydelta,y<=self.ymax+ydelta)
- x = NX.compress(testy, x)
- y = NX.compress(testy, y)
+ testx = npy.logical_and(x>=self.xmin-xdelta,x<=self.xmax+xdelta)
+ x = npy.compress(testx, x)
+ y = npy.compress(testx, y)
+ testy = npy.logical_and(y>=self.ymin-ydelta,y<=self.ymax+ydelta)
+ x = npy.compress(testy, x)
+ y = npy.compress(testy, y)
if len(x) > 1 and len(y) > 1:
# split into separate line segments if necessary.
# (not necessary for mercator or cylindrical or miller).
xd = (x[1:]-x[0:-1])**2
yd = (y[1:]-y[0:-1])**2
- dist = NX.sqrt(xd+yd)
+ dist = npy.sqrt(xd+yd)
split = dist > 500000.
- if NX.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
- ind = (NX.compress(split,squeeze(split*NX.indices(xd.shape)))+1).tolist()
+ if npy.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
+ ind = (npy.compress(split,squeeze(split*npy.indices(xd.shape)))+1).tolist()
xl = []
yl = []
iprev = 0
@@ -1564,10 +1640,10 @@
nmax = int((self.ymax-self.ymin)/dy+1)
yy = linspace(self.llcrnry,self.urcrnry,nmax)
if side == 'l':
- lons,lats = self(self.llcrnrx*NX.ones(yy.shape,NX.Float32),yy,inverse=True)
+ lons,lats = self(self.llcrnrx*npy.ones(yy.shape,npy.float32),yy,inverse=True)
lons = lons.tolist(); lats = lats.tolist()
else:
- lons,lats = self(self.urcrnrx*NX.ones(yy.shape,NX.Float32),yy,inverse=True)
+ lons,lats = self(self.urcrnrx*npy.ones(yy.shape,npy.float32),yy,inverse=True)
lons = lons.tolist(); lats = lats.tolist()
if max(lons) > 1.e20 or max(lats) > 1.e20:
raise ValueError,'inverse transformation undefined - please adjust the map projection region'
@@ -1577,10 +1653,10 @@
nmax = int((self.xmax-self.xmin)/dx+1)
xx = linspace(self.llcrnrx,self.urcrnrx,nmax)
if side == 'b':
- lons,lats = self(xx,self.llcrnry*NX.ones(xx.shape,NX.Float32),inverse=True)
+ lons,lats = self(xx,self.llcrnry*npy.ones(xx.shape,npy.float32),inverse=True)
lons = lons.tolist(); lats = lats.tolist()
else:
- lons,lats = self(xx,self.urcrnry*NX.ones(xx.shape,NX.Float32),inverse=True)
+ lons,lats = self(xx,self.urcrnry*npy.ones(xx.shape,npy.float32),inverse=True)
lons = lons.tolist(); lats = lats.tolist()
if max(lons) > 1.e20 or max(lats) > 1.e20:
raise ValueError,'inverse transformation undefined - please adjust the map projection region'
@@ -1605,7 +1681,7 @@
lonlabstr = u'-%s\N{DEGREE SIGN}'%fmt
else:
lonlabstr = u'%s\N{DEGREE SIGN}W'%fmt
- lonlab = lonlabstr%NX.fabs(lon-360)
+ lonlab = lonlabstr%npy.fabs(lon-360)
elif lon<180 and lon != 0:
if rcParams['text.usetex']:
if labelstyle=='+/-':
@@ -1721,8 +1797,8 @@
raise ValueError, 'lons and lats must be increasing!'
# check that lons in -180,180 for non-cylindrical projections.
if self.projection not in ['cyl','merc','mill']:
- lonsa = NX.array(lons)
- count = NX.sum(lonsa < -180.00001) + NX.sum(lonsa > 180.00001)
+ lonsa = npy.array(lons)
+ count = npy.sum(lonsa < -180.00001) + npy.sum(lonsa > 180.00001)
if count > 1:
raise ValueError,'grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)'
# allow for wraparound point to be outside.
@@ -1775,8 +1851,8 @@
raise ValueError, 'lons and lats must be increasing!'
# check that lons in -180,180 for non-cylindrical projections.
if self.projection not in ['cyl','merc','mill']:
- lonsa = NX.array(lons)
- count = NX.sum(lonsa < -180.00001) + NX.sum(lonsa > 180.00001)
+ lonsa = npy.array(lons)
+ count = npy.sum(lonsa < -180.00001) + npy.sum(lonsa > 180.00001)
if count > 1:
raise ValueError,'grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)'
# allow for wraparound point to be outside.
@@ -2044,8 +2120,8 @@
ax = popd(kwargs,'ax')
# make x,y masked arrays
# (masked where data is outside of projection limb)
- x = ma.masked_values(NX.where(x > 1.e20,1.e20,x), 1.e20)
- y = ma.masked_values(NX.where(y > 1.e20,1.e20,y), 1.e20)
+ x = ma.masked_values(npy.where(x > 1.e20,1.e20,x), 1.e20)
+ y = ma.masked_values(npy.where(y > 1.e20,1.e20,y), 1.e20)
# allow callers to override the hold state by passing hold=True|False
b = ax.ishold()
h = popd(kwargs, 'hold', None)
@@ -2147,10 +2223,10 @@
function to adjust the data to be consistent with the map projection
region (see examples/contour_demo.py)."""
# mask for points outside projection limb.
- xymask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))
+ xymask = npy.logical_or(npy.greater(x,1.e20),npy.greater(y,1.e20))
data = ma.asarray(data)
# combine with data mask.
- mask = NX.logical_or(ma.getmaskarray(data),xymask)
+ mask = npy.logical_or(ma.getmaskarray(data),xymask)
data = ma.masked_array(data,mask=mask)
# allow callers to override the hold state by passing hold=True|False
b = ax.ishold()
@@ -2215,10 +2291,10 @@
function to adjust the data to be consistent with the map projection
region (see examples/contour_demo.py)."""
# mask for points outside projection limb.
- xymask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))
+ xymask = npy.logical_or(npy.greater(x,1.e20),npy.greater(y,1.e20))
data = ma.asarray(data)
# combine with data mask.
- mask = NX.logical_or(ma.getmaskarray(data),xymask)
+ mask = npy.logical_or(ma.getmaskarray(data),xymask)
data = ma.masked_array(data,mask=mask)
# allow callers to override the hold state by passing hold=True|False
b = ax.ishold()
@@ -2339,9 +2415,9 @@
lsmaskf = open(os.path.join(basemap_datadir,'5minmask.bin'),'rb')
nlons = 4320; nlats = nlons/2
delta = 360./float(nlons)
- lsmask_lons = NX.arange(-180+0.5*delta,180.,delta)
- lsmask_lats = NX.arange(-90.+0.5*delta,90.,delta)
- lsmask = NX.reshape(NX.fromstring(lsmaskf.read(),NX.UInt8),(nlats,nlons))
+ lsmask_lons = npy.arange(-180+0.5*delta,180.,delta)
+ lsmask_lats = npy.arange(-90.+0.5*delta,90.,delta)
+ lsmask = npy.reshape(npy.fromstring(lsmaskf.read(),npy.uint8),(nlats,nlons))
lsmaskf.close()
# instance variable lsmask is set on first invocation,
# it contains the land-sea mask interpolated to the native
@@ -2367,17 +2443,17 @@
self.lsmask = mask
# optionally, set lakes to ocean color.
if lakes:
- mask = NX.where(self.lsmask==2,0,self.lsmask)
+ mask = npy.where(self.lsmask==2,0,self.lsmask)
else:
mask = self.lsmask
ny, nx = mask.shape
- rgba = NX.ones((ny,nx,4),NX.UInt8)
- rgba_land = NX.array(rgba_land,NX.UInt8)
- rgba_ocean = NX.array(rgba_ocean,NX.UInt8)
+ rgba = npy.ones((ny,nx,4),npy.uint8)
+ rgba_land = npy.array(rgba_land,npy.uint8)
+ rgba_ocean = npy.array(rgba_ocean,npy.uint8)
for k in range(4):
- rgba[:,:,k] = NX.where(mask,rgba_land[k],rgba_ocean[k])
+ rgba[:,:,k] = npy.where(mask,rgba_land[k],rgba_ocean[k])
# make points outside projection limb transparent.
- rgba[:,:,3] = NX.where(mask==255,0,rgba[:,:,3])
+ rgba[:,:,3] = npy.where(mask==255,0,rgba[:,:,3])
# plot mask as rgba image.
im = self.imshow(rgba,interpolation='nearest',ax=ax,**kwargs)
@@ -2455,10 +2531,10 @@
# check that xout,yout are
# within region defined by xin,yin.
if checkbounds:
- if min(NX.ravel(xout)) < min(xin) or \
- max(NX.ravel(xout)) > max(xin) or \
- min(NX.ravel(yout)) < min(yin) or \
- max(NX.ravel(yout)) > max(yin):
+ if xout.min() < xin.min() or \
+ xout.max() > xin.max() or \
+ yout.min() < yin.min() or \
+ yout.max() > yin.max():
raise ValueError, 'yout or xout outside range of yin or xin'
# compute grid coordinates of output grid.
delx = xin[1:]-xin[0:-1]
@@ -2469,9 +2545,9 @@
ycoords = (len(yin)-1)*(yout-yin[0])/(yin[-1]-yin[0])
else:
# irregular (but still rectilinear) input grid.
- xoutflat = NX.ravel(xout); youtflat = NX.ravel(yout)
- ix = (NX.searchsorted(xin,xoutflat)-1).tolist()
- iy = (NX.searchsorted(yin,youtflat)-1).tolist()
+ xoutflat = xout.flatten(); youtflat = yout.flatten()
+ ix = (npy.searchsorted(xin,xoutflat)-1).tolist()
+ iy = (npy.searchsorted(yin,youtflat)-1).tolist()
xoutflat = xoutflat.tolist(); xin = xin.tolist()
youtflat = youtflat.tolist(); yin = yin.tolist()
xcoords = []; ycoords = []
@@ -2489,33 +2565,33 @@
ycoords.append(len(yin)) # outside range on upper end
else:
ycoords.append(float(j)+(youtflat[m]-yin[j])/(yin[j+1]-yin[j]))
- xcoords = NX.reshape(xcoords,xout.shape)
- ycoords = NX.reshape(ycoords,yout.shape)
+ xcoords = npy.reshape(xcoords,xout.shape)
+ ycoords = npy.reshape(ycoords,yout.shape)
# data outside range xin,yin will be clipped to
# values on boundary.
if masked:
- xmask = NX.logical_or(NX.less(xcoords,0),NX.greater(xcoords,len(xin)-1))
- ymask = NX.logical_or(NX.less(ycoords,0),NX.greater(ycoords,len(yin)-1))
- xymask = NX.logical_or(xmask,ymask)
- xcoords = NX.clip(xcoords,0,len(xin)-1)
- ycoords = NX.clip(ycoords,0,len(yin)-1)
+ xmask = npy.logical_or(npy.less(xcoords,0),npy.greater(xcoords,len(xin)-1))
+ ymask = npy.logical_or(npy.less(ycoords,0),npy.greater(ycoords,len(yin)-1))
+ xymask = npy.logical_or(xmask,ymask)
+ xcoords = npy.clip(xcoords,0,len(xin)-1)
+ ycoords = npy.clip(ycoords,0,len(yin)-1)
# interpolate to output grid using bilinear interpolation.
if order == 1:
- xi = xcoords.astype(NX.Int32)
- yi = ycoords.astype(NX.Int32)
+ xi = xcoords.astype(npy.int32)
+ yi = ycoords.astype(npy.int32)
xip1 = xi+1
yip1 = yi+1
- xip1 = NX.clip(xip1,0,len(xin)-1)
- yip1 = NX.clip(yip1,0,len(yin)-1)
- delx = xcoords-xi.astype(NX.Float32)
- dely = ycoords-yi.astype(NX.Float32)
+ xip1 = npy.clip(xip1,0,len(xin)-1)
+ yip1 = npy.clip(yip1,0,len(yin)-1)
+ delx = xcoords-xi.astype(npy.float32)
+ dely = ycoords-yi.astype(npy.float32)
dataout = (1.-delx)*(1.-dely)*datain[yi,xi] + \
delx*dely*datain[yip1,xip1] + \
(1.-delx)*dely*datain[yip1,xi] + \
delx*(1.-dely)*datain[yi,xip1]
elif order == 0:
- xcoordsi = NX.around(xcoords).astype(NX.Int32)
- ycoordsi = NX.around(ycoords).astype(NX.Int32)
+ xcoordsi = npy.around(xcoords).astype(npy.int32)
+ ycoordsi = npy.around(ycoords).astype(npy.int32)
dataout = datain[ycoordsi,xcoordsi]
else:
raise ValueError,'order keyword must be 0 or 1'
@@ -2524,7 +2600,7 @@
newmask = ma.mask_or(ma.getmask(dataout), xymask)
dataout = ma.masked_array(dataout,mask=newmask)
elif masked and is_scalar(masked):
- dataout = NX.where(xymask,masked,dataout)
+ dataout = npy.where(xymask,masked,dataout)
return dataout
def shiftgrid(lon0,datain,lonsin,start=True):
@@ -2542,13 +2618,13 @@
returns dataout,lonsout (data and longitudes on shifted grid).
"""
- if NX.fabs(lonsin[-1]-lonsin[0]-360.) > 1.e-4:
+ if npy.fabs(lonsin[-1]-lonsin[0]-360.) > 1.e-4:
raise ValueError, 'cyclic point not included'
if lon0 < lonsin[0] or lon0 > lonsin[-1]:
raise ValueError, 'lon0 outside of range of lonsin'
- i0 = NX.argmin(NX.fabs(lonsin-lon0))
- dataout = NX.zeros(datain.shape,datain.dtype)
- lonsout = NX.zeros(lonsin.shape,lonsin.dtype)
+ i0 = npy.argmin(npy.fabs(lonsin-lon0))
+ dataout = npy.zeros(datain.shape,datain.dtype)
+ lonsout = npy.zeros(lonsin.shape,lonsin.dtype)
if start:
lonsout[0:len(lonsin)-i0] = lonsin[i0:]
else:
@@ -2572,13 +2648,13 @@
if hasattr(arrin,'mask'):
arrout = ma.zeros((nlats,nlons+1),arrin.dtype)
else:
- arrout = NX.zeros((nlats,nlons+1),arrin.dtype)
+ arrout = npy.zeros((nlats,nlons+1),arrin.dtype)
arrout[:,0:nlons] = arrin[:,:]
arrout[:,nlons] = arrin[:,0]
if hasattr(lonsin,'mask'):
lonsout = ma.zeros(nlons+1,lonsin.dtype)
else:
- lonsout = NX.zeros(nlons+1,lonsin.dtype)
+ lonsout = npy.zeros(nlons+1,lonsin.dtype)
lonsout[0:nlons] = lonsin[:]
lonsout[nlons] = lonsin[-1] + lonsin[1]-lonsin[0]
return arrout,lonsout
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-03 12:22:45
|
Revision: 4099
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4099&view=rev
Author: jswhit
Date: 2007-11-03 05:22:43 -0700 (Sat, 03 Nov 2007)
Log Message:
-----------
only read in state, country and river geometries when draw method invoked.
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-02 23:21:06 UTC (rev 4098)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-03 12:22:43 UTC (rev 4099)
@@ -661,9 +661,6 @@
# set defaults for area_thresh.
self.resolution = resolution
- # if no boundary data needed, we are done.
- if self.resolution is None:
- return
if area_thresh is None:
if resolution == 'c':
area_thresh = 10000.
@@ -680,18 +677,14 @@
self._boundarypolyll, self._boundarypolyxy = self._getmapboundary()
# read in coastline polygons, only keeping those that
# intersect map boundary polygon.
- self.coastsegs, self.coastpolygontypes = self._readboundarydata('gshhs')
- # same for countries, states, rivers.
- self.cntrysegs, types = self._readboundarydata('countries')
- self.statesegs, types = self._readboundarydata('states')
- self.riversegs, types = self._readboundarydata('rivers')
- # for coastlines, reformat for use in
- # matplotlib.patches.Polygon.
- self.coastpolygons = []
- for xy in self.coastsegs:
- x = [x1 for x1,x2 in xy]
- y = [x2 for x1,x2 in xy]
- self.coastpolygons.append((x,y))
+ if self.resolution is not None:
+ self.coastsegs, self.coastpolygontypes = self._readboundarydata('gshhs')
+ # reformat for use in matplotlib.patches.Polygon.
+ self.coastpolygons = []
+ for xy in self.coastsegs:
+ x = [x1 for x1,x2 in xy]
+ y = [x2 for x1,x2 in xy]
+ self.coastpolygons.append((x,y))
def __call__(self,x,y,inverse=False):
"""
@@ -791,22 +784,25 @@
else:
continue
if poly.intersects(self._boundarypolyll):
- poly = poly.intersection(self._boundarypolyll)
- if hasattr(poly,'geoms'):
- geoms = poly.geoms
- else:
- geoms = [poly]
- for psub in geoms:
- if name == 'gshhs':
- b = npy.asarray(psub.boundary)
+ try:
+ poly = poly.intersection(self._boundarypolyll)
+ if hasattr(poly,'geoms'):
+ geoms = poly.geoms
else:
- b = npy.asarray(psub.coords)
- blons = b[:,0]; blats = b[:,1]
- bx, by = self(blons, blats)
- #if (bx > 1.20).any() or (by > 1.e20).any():
- # continue
- polygons.append(zip(bx,by))
- polygon_types.append(type)
+ geoms = [poly]
+ for psub in geoms:
+ if name == 'gshhs':
+ b = npy.asarray(psub.boundary)
+ else:
+ b = npy.asarray(psub.coords)
+ blons = b[:,0]; blats = b[:,1]
+ bx, by = self(blons, blats)
+ #if (bx > 1.20).any() or (by > 1.e20).any():
+ # continue
+ polygons.append(zip(bx,by))
+ polygon_types.append(type)
+ except:
+ pass
# if map boundary polygon is not valid in lat/lon
# coordinates, compute intersection between map
# projection region and boundary geometries in map
@@ -860,7 +856,6 @@
b = npy.asarray(psub.coords)
polygons.append(zip(b[:,0],b[:,1]))
polygon_types.append(type)
- print name,len(polygons)
return polygons, polygon_types
@@ -1119,6 +1114,10 @@
"""
if self.resolution is None:
raise AttributeError, 'there are no boundary datasets associated with this Basemap instance'
+ # read in country line segments, only keeping those that
+ # intersect map boundary polygon.
+ if not hasattr(self,'cntrysegs'):
+ self.cntrysegs, types = self._readboundarydata('countries')
# get current axes instance (if none specified).
if ax is None and self.ax is None:
try:
@@ -1150,6 +1149,10 @@
"""
if self.resolution is None:
raise AttributeError, 'there are no boundary datasets associated with this Basemap instance'
+ # read in state line segments, only keeping those that
+ # intersect map boundary polygon.
+ if not hasattr(self,'statesegs'):
+ self.statesegs, types = self._readboundarydata('states')
# get current axes instance (if none specified).
if ax is None and self.ax is None:
try:
@@ -1181,6 +1184,10 @@
"""
if self.resolution is None:
raise AttributeError, 'there are no boundary datasets associated with this Basemap instance'
+ # read in river line segments, only keeping those that
+ # intersect map boundary polygon.
+ if not hasattr(self,'riversegs'):
+ self.riversegs, types = self._readboundarydata('rivers')
# get current axes instance (if none specified).
if ax is None and self.ax is None:
try:
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-04 15:57:22
|
Revision: 4103
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4103&view=rev
Author: jswhit
Date: 2007-11-04 07:57:19 -0800 (Sun, 04 Nov 2007)
Log Message:
-----------
remove debug code
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-04 06:09:22 UTC (rev 4102)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-04 15:57:19 UTC (rev 4103)
@@ -777,10 +777,6 @@
lons.append(lonend)
lats.append(-90.)
poly = PolygonShape(zip(lons,lats))
- #b = npy.asarray(poly.boundary)
- #import pylab
- #pylab.fill(b[:,0],b[:,1],'b')
- #pylab.show()
else:
continue
if poly.intersects(self._boundarypolyll):
@@ -797,8 +793,6 @@
b = npy.asarray(psub.coords)
blons = b[:,0]; blats = b[:,1]
bx, by = self(blons, blats)
- #if (bx > 1.20).any() or (by > 1.e20).any():
- # continue
polygons.append(zip(bx,by))
polygon_types.append(type)
except:
@@ -814,6 +808,8 @@
b = npy.asarray(poly.coords)
bx, by = self(b[:,0], b[:,1])
mask = npy.logical_and(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 name == 'gshhs':
poly = PolygonShape(zip(bx,by))
@@ -823,23 +819,6 @@
bx = npy.compress(mask, bx)
by = npy.compress(mask, by)
poly = LineShape(zip(bx,by))
- if 0:
- import pylab
- print poly
- pp = self._boundarypolyxy.intersection(poly)
- print name, pp
- a = npy.asarray(self._boundarypolyxy.boundary)
- pylab.plot(a[:,0],a[:,1],'b')
- pylab.plot(bx,by,'r')
- if hasattr(pp,'geoms'):
- px = pp.geoms
- else:
- px = [pp]
- for p in px:
- c = npy.asarray(p.boundary)
- pylab.fill(c[:,0],c[:,1],'g')
- pylab.show()
- raise SystemExit(0)
if self._boundarypolyxy.intersects(poly):
try:
poly = self._boundarypolyxy.intersection(poly)
@@ -858,7 +837,6 @@
polygon_types.append(type)
return polygons, polygon_types
-
def _getmapboundary(self):
"""
define map boundary polygon (in lat/lon coordinates)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-04 22:14:43
|
Revision: 4104
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4104&view=rev
Author: jswhit
Date: 2007-11-04 14:14:41 -0800 (Sun, 04 Nov 2007)
Log Message:
-----------
fix Antarctica for moll,robin and sinu (moll still broken)
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-04 15:57:19 UTC (rev 4103)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-04 22:14:41 UTC (rev 4104)
@@ -758,7 +758,7 @@
# regions and high-resolution boundary geometries).
if not containsPole:
# close Antarctica for cylindrical projections.
- if name == 'gshhs' and self.projection in ['cyl','merc','mill']:
+ if name == 'gshhs' and self.projection in ['cyl','merc','mill','moll','robin','sinu']:
b = npy.asarray(poly.boundary)
lons = b[:,0]
lats = b[:,1]
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-05 01:36:18
|
Revision: 4105
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4105&view=rev
Author: jswhit
Date: 2007-11-04 17:36:15 -0800 (Sun, 04 Nov 2007)
Log Message:
-----------
fixed mollweide
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-04 22:14:41 UTC (rev 4104)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-05 01:36:15 UTC (rev 4105)
@@ -862,26 +862,28 @@
# quasi-elliptical region.
x = []; y = []
# left side
- lats = linspace(-89.9,89.9,ny).tolist()
- lons = len(lats)*[self.projparams['lon_0']-179.9]
- x,y = self(lons,lats)
+ lats1 = linspace(-89.9,89.9,ny).tolist()
+ lons1 = len(lats1)*[self.projparams['lon_0']-179.9]
+ x,y = self(lons1,lats1)
# top.
- lons = linspace(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179,nx).tolist()
- lats = len(lons)*[89.9]
- xx,yy = self(lons,lats)
+ lons2 = linspace(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179,nx).tolist()
+ lats2 = len(lons2)*[89.9]
+ xx,yy = self(lons2,lats2)
x = x+xx; y = y+yy
# right side
- lats = linspace(89.9,-89.9,ny).tolist()
- lons = len(lats)*[self.projparams['lon_0']+179.9]
- xx,yy = self(lons,lats)
+ lats3 = linspace(89.9,-89.9,ny).tolist()
+ lons3 = len(lats3)*[self.projparams['lon_0']+179.9]
+ xx,yy = self(lons3,lats3)
x = x+xx; y = y+yy
# bottom.
- lons = linspace(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180).tolist()
- lats = len(lons)*[-89.9]
- xx,yy = self(lons,lats)
+ lons4 = linspace(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180).tolist()
+ lats4 = len(lons4)*[-89.9]
+ xx,yy = self(lons4,lats4)
x = x+xx; y = y+yy
x = npy.array(x,npy.float64)
y = npy.array(y,npy.float64)
+ lons = lons1+lons2+lons3+lons4
+ lats = lats1+lats2+lats3+lats4
boundaryxy = PolygonShape(zip(x,y))
else: # all other projections are rectangular.
# left side (x = xmin, ymin <= y <= ymax)
@@ -908,21 +910,21 @@
lons = [self.llcrnrlon, self.llcrnrlon, self.urcrnrlon, self.urcrnrlon]
lats = [self.llcrnrlat, self.urcrnrlat, self.urcrnrlat, self.llcrnrlat]
else:
- lons, lats = self(x,y,inverse=True)
- # fix lons so there are no jumps.
- n = 1
- lonprev = lons[0]
- for lon,lat in zip(lons[1:],lats[1:]):
- if npy.abs(lon-lonprev) > 90.:
- if lonprev < 0:
- lon = lon - 360.
- else:
- lon = lon + 360
- lons[n] = lon
- lonprev = lon
- n = n + 1
+ if self.projection not in ['moll','robin','sinu']:
+ lons, lats = self(x,y,inverse=True)
+ # fix lons so there are no jumps.
+ n = 1
+ lonprev = lons[0]
+ for lon,lat in zip(lons[1:],lats[1:]):
+ if npy.abs(lon-lonprev) > 90.:
+ if lonprev < 0:
+ lon = lon - 360.
+ else:
+ lon = lon + 360
+ lons[n] = lon
+ lonprev = lon
+ n = n + 1
boundaryll = PolygonShape(zip(lons,lats))
- #print 'map projection region',boundaryll.geom_type,boundaryll.is_valid
return PolygonShape(zip(lons,lats)), boundaryxy
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-05 16:11:30
|
Revision: 4108
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4108&view=rev
Author: jswhit
Date: 2007-11-05 08:11:12 -0800 (Mon, 05 Nov 2007)
Log Message:
-----------
fixes for non-fulldisk geos and ortho, cleanups in _getmapboundary
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-05 15:45:00 UTC (rev 4107)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-05 16:11:12 UTC (rev 4108)
@@ -734,6 +734,10 @@
hasNP = self._boundarypolyxy.contains(NPole)
hasSP = self._boundarypolyxy.contains(SPole)
containsPole = hasNP or hasSP
+ # make sure geostationary projection has containsPole=False
+ # even if only a limited area is requested (since entire
+ # disk will be used to subset boundary geometries).
+ if self.projection == 'geos': containsPole = False
# iterate over boundary geometries.
for line in bdatmetafile:
linesplit = line.split()
@@ -843,47 +847,55 @@
"""
nx = 100
ny = 100
- if self.projection == 'ortho' and self._fulldisk:
+ maptran = self
+ if self.projection in ['ortho','geos']:
# circular region.
thetas = linspace(0.,2.*npy.pi,2*nx*ny)[:-1]
- radius = self.rmajor
- x = radius*npy.cos(thetas) + 0.5*self.xmax
- y = radius*npy.sin(thetas) + 0.5*self.ymax
+ if self.projection == 'ortho':
+ rminor = self.rmajor
+ rmajor = self.rmajor
+ else:
+ rminor = self._height
+ rmajor = self._width
+ x = rmajor*npy.cos(thetas) + rmajor
+ y = rminor*npy.sin(thetas) + rminor
boundaryxy = PolygonShape(zip(x,y))
- elif self.projection == 'geos' and self._fulldisk:
- # elliptical region
- thetas = linspace(0.,2.*npy.pi,2*nx*ny)[:-1]
- rminor = self._height
- rmajor = self._width
- x = rmajor*npy.cos(thetas) + 0.5*self.xmax
- y = rminor*npy.sin(thetas) + 0.5*self.ymax
- boundaryxy = PolygonShape(zip(x,y))
+ # compute proj instance for full disk, if necessary.
+ if not self._fulldisk:
+ projparms = self.projparams
+ del projparms['x_0']
+ del projparms['y_0']
+ if self.projection == 'ortho':
+ llcrnrx = -self.rmajor
+ llcrnry = -self.rmajor
+ urcrnrx = -llcrnrx
+ urcrnry = -llcrnry
+ else:
+ llcrnrx = -self._width
+ llcrnry = -self._height
+ urcrnrx = -llcrnrx
+ urcrnry = -llcrnry
+ projparms['x_0']=-llcrnrx
+ projparms['y_0']=-llcrnry
+ maptran = pyproj.Proj(projparms)
elif self.projection in ['moll','robin','sinu']:
# quasi-elliptical region.
- x = []; y = []
+ lon_0 = self.projparams['lon_0']
# left side
lats1 = linspace(-89.9,89.9,ny).tolist()
- lons1 = len(lats1)*[self.projparams['lon_0']-179.9]
- x,y = self(lons1,lats1)
+ lons1 = len(lats1)*[lon_0-179.9]
# top.
- lons2 = linspace(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179,nx).tolist()
+ lons2 = linspace(lon_0-179.9,lon_0+179.9,nx).tolist()
lats2 = len(lons2)*[89.9]
- xx,yy = self(lons2,lats2)
- x = x+xx; y = y+yy
# right side
lats3 = linspace(89.9,-89.9,ny).tolist()
- lons3 = len(lats3)*[self.projparams['lon_0']+179.9]
- xx,yy = self(lons3,lats3)
- x = x+xx; y = y+yy
+ lons3 = len(lats3)*[lon_0+179.9]
# bottom.
- lons4 = linspace(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180).tolist()
+ lons4 = linspace(lon_0+179.9,lon_0-179.9,nx).tolist()
lats4 = len(lons4)*[-89.9]
- xx,yy = self(lons4,lats4)
- x = x+xx; y = y+yy
- x = npy.array(x,npy.float64)
- y = npy.array(y,npy.float64)
- lons = lons1+lons2+lons3+lons4
- lats = lats1+lats2+lats3+lats4
+ lons = npy.array(lons1+lons2+lons3+lons4,npy.float64)
+ lats = npy.array(lats1+lats2+lats3+lats4,npy.float64)
+ x, y = maptran(lons,lats)
boundaryxy = PolygonShape(zip(x,y))
else: # all other projections are rectangular.
# left side (x = xmin, ymin <= y <= ymax)
@@ -911,7 +923,7 @@
lats = [self.llcrnrlat, self.urcrnrlat, self.urcrnrlat, self.llcrnrlat]
else:
if self.projection not in ['moll','robin','sinu']:
- lons, lats = self(x,y,inverse=True)
+ lons, lats = maptran(x,y,inverse=True)
# fix lons so there are no jumps.
n = 1
lonprev = lons[0]
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-05 16:50:25
|
Revision: 4109
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4109&view=rev
Author: jswhit
Date: 2007-11-05 08:50:11 -0800 (Mon, 05 Nov 2007)
Log Message:
-----------
fix for coastline segments that jump across plot.
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-05 16:11:12 UTC (rev 4108)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-05 16:50:11 UTC (rev 4109)
@@ -682,9 +682,27 @@
# reformat for use in matplotlib.patches.Polygon.
self.coastpolygons = []
for xy in self.coastsegs:
- x = [x1 for x1,x2 in xy]
- y = [x2 for x1,x2 in xy]
+ x, y = zip(*xy)
self.coastpolygons.append((x,y))
+ # split coastline segments that jump across entire plot.
+ coastsegs = []
+ for seg in self.coastsegs:
+ x, y = zip(*seg)
+ x = npy.array(x,npy.float64); y = npy.array(y,npy.float64)
+ xd = (x[1:]-x[0:-1])**2
+ yd = (y[1:]-y[0:-1])**2
+ dist = npy.sqrt(xd+yd)
+ split = dist > 5000000.
+ if npy.sum(split) and self.projection not in ['merc','cyl','mill']:
+ ind = (npy.compress(split,squeeze(split*npy.indices(xd.shape)))+1).tolist()
+ iprev = 0
+ ind.append(len(xd))
+ for i in ind:
+ coastsegs.append(zip(x[iprev:i],y[iprev:i]))
+ iprev = i
+ else:
+ coastsegs.append(seg)
+ self.coastsegs = coastsegs
def __call__(self,x,y,inverse=False):
"""
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-05 19:47:41
|
Revision: 4111
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4111&view=rev
Author: jswhit
Date: 2007-11-05 11:47:27 -0800 (Mon, 05 Nov 2007)
Log Message:
-----------
support for orthographic
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-05 17:30:08 UTC (rev 4110)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-05 19:47:27 UTC (rev 4111)
@@ -661,7 +661,7 @@
# set defaults for area_thresh.
self.resolution = resolution
- if area_thresh is None:
+ if area_thresh is None and resolution is not None:
if resolution == 'c':
area_thresh = 10000.
elif resolution == 'l':
@@ -681,13 +681,11 @@
self.coastsegs, self.coastpolygontypes = self._readboundarydata('gshhs')
# reformat for use in matplotlib.patches.Polygon.
self.coastpolygons = []
- for xy in self.coastsegs:
- x, y = zip(*xy)
- self.coastpolygons.append((x,y))
- # split coastline segments that jump across entire plot.
+ # also, split coastline segments that jump across entire plot.
coastsegs = []
for seg in self.coastsegs:
x, y = zip(*seg)
+ self.coastpolygons.append((x,y))
x = npy.array(x,npy.float64); y = npy.array(y,npy.float64)
xd = (x[1:]-x[0:-1])**2
yd = (y[1:]-y[0:-1])**2
@@ -749,13 +747,29 @@
# see if map projection region polygon contains a pole.
NPole = PointShape(self(0.,90.))
SPole = PointShape(self(0.,-90.))
- hasNP = self._boundarypolyxy.contains(NPole)
- hasSP = self._boundarypolyxy.contains(SPole)
+ boundarypolyxy = self._boundarypolyxy
+ hasNP = boundarypolyxy.contains(NPole)
+ hasSP = boundarypolyxy.contains(SPole)
containsPole = hasNP or hasSP
# make sure geostationary projection has containsPole=False
# even if only a limited area is requested (since entire
# disk will be used to subset boundary geometries).
if self.projection == 'geos': containsPole = False
+ # make sure orthographic projection has containsPole=True
+ # we will compute the intersections in stereographic
+ # coordinates, then transform to orthographic.
+ if self.projection == 'ortho':
+ containsPole = True
+ lon_0=self.projparams['lon_0']
+ lat_0=self.projparams['lat_0']
+ # FIXME: won't work for points exactly on equator??
+ if npy.abs(lat_0) < 1.e-4: lat_0 = 1.e04
+ maptran = pyproj.Proj(proj='stere',lon_0=lon_0,lat_0=lat_0)
+ # boundary polygon for orthographic projection
+ # in stereographic coorindates.
+ b = npy.asarray(self._boundarypolyll.boundary)
+ blons = b[:,0]; blats = b[:,1]
+ boundarypolyxy = PolygonShape(zip(*maptran(blons,blats)))
# iterate over boundary geometries.
for line in bdatmetafile:
linesplit = line.split()
@@ -828,7 +842,13 @@
b = npy.asarray(poly.boundary)
else:
b = npy.asarray(poly.coords)
- bx, by = self(b[:,0], b[:,1])
+ blons = b[:,0]; blats = b[:,1]
+ # special case for ortho, compute polygon
+ # coordinates in stereographic coords.
+ if self.projection == 'ortho':
+ bx, by = maptran(blons, blats)
+ else:
+ bx, by = self(blons, blats)
mask = npy.logical_and(bx<1.e20,by<1.e20)
# if less than two points are valid in
# map proj coords, skip this geometry.
@@ -841,9 +861,9 @@
bx = npy.compress(mask, bx)
by = npy.compress(mask, by)
poly = LineShape(zip(bx,by))
- if self._boundarypolyxy.intersects(poly):
+ if boundarypolyxy.intersects(poly):
try:
- poly = self._boundarypolyxy.intersection(poly)
+ poly = boundarypolyxy.intersection(poly)
except:
continue
if hasattr(poly,'geoms'):
@@ -855,6 +875,19 @@
b = npy.asarray(psub.boundary)
else:
b = npy.asarray(psub.coords)
+ # if projection == 'ortho',
+ # transform polygon from stereographic
+ # to orthographic coordinates.
+ if self.projection == 'ortho':
+ # if coastline polygon covers more than 99%
+ # of map region, it's probably bogus, so
+ # skip it.
+ if name == 'gshhs' and \
+ psub.area/boundarypolyxy.area > 0.99: continue
+ bx = b[:,0]; by = b[:,1]
+ blons, blats = maptran(bx, by, inverse=True)
+ bx, by = self(blons, blats)
+ b[:,0] = bx; b[:,1] = by
polygons.append(zip(b[:,0],b[:,1]))
polygon_types.append(type)
return polygons, polygon_types
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-05 20:53:51
|
Revision: 4116
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4116&view=rev
Author: jswhit
Date: 2007-11-05 12:53:49 -0800 (Mon, 05 Nov 2007)
Log Message:
-----------
some optimizations for orthographic projection
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-05 20:42:23 UTC (rev 4115)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-05 20:53:49 UTC (rev 4116)
@@ -504,6 +504,9 @@
self._fulldisk = False
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
+ # FIXME: won't work for points exactly on equator??
+ if npy.abs(lat_0) < 1.e-2: lat_0 = 1.e-2
+ projparams['lat_0'] = lat_0
elif projection == 'geos':
if lon_0 is None and satellite_height is None:
raise ValueError, 'must specify lon_0 and satellite_height for Geostationary basemap'
@@ -758,12 +761,10 @@
# make sure orthographic projection has containsPole=True
# we will compute the intersections in stereographic
# coordinates, then transform to orthographic.
- if self.projection == 'ortho':
+ if self.projection == 'ortho' and name == 'gshhs':
containsPole = True
lon_0=self.projparams['lon_0']
lat_0=self.projparams['lat_0']
- # FIXME: won't work for points exactly on equator??
- if npy.abs(lat_0) < 1.e-4: lat_0 = 1.e04
maptran = pyproj.Proj(proj='stere',lon_0=lon_0,lat_0=lat_0)
# boundary polygon for orthographic projection
# in stereographic coorindates.
@@ -845,7 +846,7 @@
blons = b[:,0]; blats = b[:,1]
# special case for ortho, compute polygon
# coordinates in stereographic coords.
- if self.projection == 'ortho':
+ if name == 'gshhs' and self.projection == 'ortho':
bx, by = maptran(blons, blats)
else:
bx, by = self(blons, blats)
@@ -860,6 +861,13 @@
# in this map projection.
bx = npy.compress(mask, bx)
by = npy.compress(mask, by)
+ # for orthographic projection, all points
+ # outside map projection region are eliminated
+ # with the above step, so we're done.
+ if self.projection == 'ortho':
+ polygons.append(zip(bx,by))
+ polygon_types.append(type)
+ continue
poly = LineShape(zip(bx,by))
if boundarypolyxy.intersects(poly):
try:
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-05 21:02:36
|
Revision: 4117
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4117&view=rev
Author: jswhit
Date: 2007-11-05 13:02:33 -0800 (Mon, 05 Nov 2007)
Log Message:
-----------
cleanups in drawmapboundary
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-05 20:53:49 UTC (rev 4116)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-05 21:02:33 UTC (rev 4117)
@@ -19,7 +19,7 @@
from shapely.geometry import Polygon as PolygonShape
from shapely.geometry import LineString as LineShape
from shapely.geometry import Point as PointShape
-from shapely import wkb, wkt
+from shapely import wkb
# basemap data files now installed in lib/matplotlib/toolkits/basemap/data
basemap_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
@@ -1013,9 +1013,6 @@
ax = pylab.gca()
elif ax is None and self.ax is not None:
ax = self.ax
- x = []
- y = []
- dtheta = 0.01
if self.projection == 'ortho' and self._fulldisk: # circular region.
# define a circle patch, add it to axes instance.
circle = Circle((self.rmajor,self.rmajor),self.rmajor)
@@ -1033,25 +1030,24 @@
ellps.set_linewidth(linewidth)
ellps.set_clip_on(False)
elif self.projection in ['moll','robin','sinu']: # elliptical region.
+ nx = 100; ny = 100
+ # quasi-elliptical region.
+ lon_0 = self.projparams['lon_0']
# left side
- lats = npy.arange(-89.9,89.9+dtheta,dtheta).tolist()
- lons = len(lats)*[self.projparams['lon_0']-179.9]
- x,y = self(lons,lats)
+ lats1 = linspace(-89.9,89.9,ny).tolist()
+ lons1 = len(lats1)*[lon_0-179.9]
# top.
- lons = npy.arange(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179+dtheta,dtheta).tolist()
- lats = len(lons)*[89.9]
- xx,yy = self(lons,lats)
- x = x+xx; y = y+yy
+ lons2 = linspace(lon_0-179.9,lon_0+179.9,nx).tolist()
+ lats2 = len(lons2)*[89.9]
# right side
- lats = npy.arange(89.9,-89.9-dtheta,-dtheta).tolist()
- lons = len(lats)*[self.projparams['lon_0']+179.9]
- xx,yy = self(lons,lats)
- x = x+xx; y = y+yy
+ lats3 = linspace(89.9,-89.9,ny).tolist()
+ lons3 = len(lats3)*[lon_0+179.9]
# bottom.
- lons = npy.arange(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180-dtheta,-dtheta).tolist()
- lats = len(lons)*[-89.9]
- xx,yy = self(lons,lats)
- x = x+xx; y = y+yy
+ lons4 = linspace(lon_0+179.9,lon_0-179.9,nx).tolist()
+ lats4 = len(lons4)*[-89.9]
+ lons = npy.array(lons1+lons2+lons3+lons4,npy.float64)
+ lats = npy.array(lats1+lats2+lats3+lats4,npy.float64)
+ x, y = self(lons,lats)
xy = zip(x,y)
poly = Polygon(xy,edgecolor=color,linewidth=linewidth)
ax.add_patch(poly)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-05 21:07:35
|
Revision: 4118
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4118&view=rev
Author: jswhit
Date: 2007-11-05 13:07:32 -0800 (Mon, 05 Nov 2007)
Log Message:
-----------
remove last vestiges of numerix
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-05 21:02:33 UTC (rev 4117)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-05 21:07:32 UTC (rev 4118)
@@ -10,10 +10,8 @@
from matplotlib.lines import Line2D
import pyproj, sys, os, math, dbflib
from proj import Proj
-from matplotlib.numerix import ma
import numpy as npy
-from numpy import linspace
-from matplotlib.numerix.mlab import squeeze
+from numpy import linspace, squeeze, ma
from matplotlib.cbook import popd, is_scalar
from shapelib import ShapeFile
from shapely.geometry import Polygon as PolygonShape
@@ -2572,7 +2570,7 @@
points in datain are masked. To avoid this, do the interpolation in two
passes, first with order=1 (producing dataout1), then with order=0
(producing dataout2). Then replace all the masked values in dataout1
- with the corresponding elements in dataout2 (using numerix.where).
+ with the corresponding elements in dataout2 (using numpy.where).
This effectively uses nearest neighbor interpolation if any of the
four surrounding points in datain are masked, and bilinear interpolation
otherwise.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-06 18:23:15
|
Revision: 4125
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4125&view=rev
Author: jswhit
Date: 2007-11-06 10:23:13 -0800 (Tue, 06 Nov 2007)
Log Message:
-----------
docstring cleanup, remove extra Shapely Polygon call.
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 15:50:51 UTC (rev 4124)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-06 18:23:13 UTC (rev 4125)
@@ -792,7 +792,8 @@
# coordinates (this saves time, especially for small map
# regions and high-resolution boundary geometries).
if not containsPole:
- # close Antarctica for cylindrical projections.
+ # close Antarctica for cylindrical and
+ # pseudo-cylindrical projections.
if name == 'gshhs' and self.projection in ['cyl','merc','mill','moll','robin','sinu']:
b = npy.asarray(poly.boundary)
lons = b[:,0]
@@ -900,7 +901,7 @@
def _getmapboundary(self):
"""
- define map boundary polygon (in lat/lon coordinates)
+ create map boundary polygon (in lat/lon and x/y coordinates)
"""
nx = 100
ny = 100
@@ -994,7 +995,7 @@
lonprev = lon
n = n + 1
boundaryll = PolygonShape(zip(lons,lats))
- return PolygonShape(zip(lons,lats)), boundaryxy
+ return boundaryll, boundaryxy
def drawmapboundary(self,color='k',linewidth=1.0,ax=None):
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-06 18:50:55
|
Revision: 4127
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4127&view=rev
Author: jswhit
Date: 2007-11-06 10:50:53 -0800 (Tue, 06 Nov 2007)
Log Message:
-----------
refine check for bogus polygons in fulldisk ortho projection.
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 18:32:30 UTC (rev 4126)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-06 18:50:53 UTC (rev 4127)
@@ -887,10 +887,12 @@
# to orthographic coordinates.
if self.projection == 'ortho':
# if coastline polygon covers more than 99%
- # of map region, it's probably bogus, so
- # skip it.
- if name == 'gshhs' and \
- psub.area/boundarypolyxy.area > 0.99: continue
+ # of map region for fulldisk projection,
+ # it's probably bogus, so skip it.
+ areafrac = psub.area/boundarypolyxy.area
+ if name == 'gshhs' and\
+ self._fulldisk and\
+ areafrac > 0.99: continue
bx = b[:,0]; by = b[:,1]
blons, blats = maptran(bx, by, inverse=True)
bx, by = self(blons, blats)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-06 19:22:39
|
Revision: 4129
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4129&view=rev
Author: jswhit
Date: 2007-11-06 11:22:37 -0800 (Tue, 06 Nov 2007)
Log Message:
-----------
more comments in _readboundar method
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 18:56:21 UTC (rev 4128)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-06 19:22:37 UTC (rev 4129)
@@ -815,36 +815,53 @@
poly = PolygonShape(zip(lons,lats))
else:
continue
+ # if polygon instersects map projection
+ # region, process it.
if poly.intersects(self._boundarypolyll):
try:
poly = poly.intersection(self._boundarypolyll)
- if hasattr(poly,'geoms'):
- geoms = poly.geoms
+ except:
+ continue
+ # create iterable object with geometries
+ # that intersect map region.
+ if hasattr(poly,'geoms'):
+ geoms = poly.geoms
+ else:
+ geoms = [poly]
+ # iterate over geometries in intersection.
+ for psub in geoms:
+ # only coastlines are polygons,
+ # which have a 'boundary' attribute.
+ # otherwise, use 'coords' attribute
+ # to extract coordinates.
+ if name == 'gshhs':
+ b = npy.asarray(psub.boundary)
else:
- geoms = [poly]
- for psub in geoms:
- if name == 'gshhs':
- b = npy.asarray(psub.boundary)
- else:
- b = npy.asarray(psub.coords)
- blons = b[:,0]; blats = b[:,1]
- bx, by = self(blons, blats)
- polygons.append(zip(bx,by))
- polygon_types.append(type)
- except:
- pass
+ b = npy.asarray(psub.coords)
+ blons = b[:,0]; blats = b[:,1]
+ # transformation from lat/lon to
+ # map projection coordinates.
+ bx, by = self(blons, blats)
+ polygons.append(zip(bx,by))
+ polygon_types.append(type)
# if map boundary polygon is not valid in lat/lon
# coordinates, compute intersection between map
# projection region and boundary geometries in map
# projection coordinates.
else:
+ # only coastlines are polygons,
+ # which have a 'boundary' attribute.
+ # otherwise, use 'coords' attribute
+ # to extract coordinates.
if name == 'gshhs':
b = npy.asarray(poly.boundary)
else:
b = npy.asarray(poly.coords)
blons = b[:,0]; blats = b[:,1]
- # special case for ortho, compute polygon
- # coordinates in stereographic coords.
+ # transform coordinates from lat/lon
+ # to map projection coordinates.
+ # special case for ortho, compute coastline polygon
+ # vertices in stereographic coords.
if name == 'gshhs' and self.projection == 'ortho':
bx, by = maptran(blons, blats)
else:
@@ -854,9 +871,12 @@
# map proj coords, skip this geometry.
if sum(mask) <= 1: continue
if name == 'gshhs':
+ # create a polygon object for coastline
+ # geometry.
poly = PolygonShape(zip(bx,by))
else:
- # remove parts of geometry that are undefined
+ # 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)
@@ -867,16 +887,22 @@
polygons.append(zip(bx,by))
polygon_types.append(type)
continue
+ # create a Line object for other geometries.
poly = LineShape(zip(bx,by))
+ # if polygon instersects map projection
+ # region, process it.
if boundarypolyxy.intersects(poly):
try:
poly = boundarypolyxy.intersection(poly)
except:
continue
+ # create iterable object with geometries
+ # that intersect map region.
if hasattr(poly,'geoms'):
geoms = poly.geoms
else:
geoms = [poly]
+ # iterate over geometries in intersection.
for psub in geoms:
if name == 'gshhs':
b = npy.asarray(psub.boundary)
@@ -894,7 +920,11 @@
self._fulldisk and\
areafrac > 0.99: continue
bx = b[:,0]; by = b[:,1]
+ # inverse transform from stereographic
+ # to lat/lon.
blons, blats = maptran(bx, by, inverse=True)
+ # forward transform from lat/lon to
+ # orthographic.
bx, by = self(blons, blats)
b[:,0] = bx; b[:,1] = by
polygons.append(zip(b[:,0],b[:,1]))
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-06 23:03:23
|
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.
|
|
From: <js...@us...> - 2007-11-07 12:55:41
|
Revision: 4138
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4138&view=rev
Author: jswhit
Date: 2007-11-07 04:55:37 -0800 (Wed, 07 Nov 2007)
Log Message:
-----------
some projections should never cross pole.
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 23:03:21 UTC (rev 4137)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-07 12:55:37 UTC (rev 4138)
@@ -749,9 +749,13 @@
NPole = PointShape(self(0.,90.))
SPole = PointShape(self(0.,-90.))
boundarypolyxy = self._boundarypolyxy
+ boundarypolyll = self._boundarypolyll
hasNP = boundarypolyxy.contains(NPole)
hasSP = boundarypolyxy.contains(SPole)
containsPole = hasNP or hasSP
+ # these projections cannot cross pole.
+ if containsPole and self.projection in ['tmerc','cass','omerc','merc']:
+ raise ValueError('%s projection cannot cross pole'%(self.projection))
# make sure geostationary projection has containsPole=False
# even if only a limited area is requested (since entire
# disk will be used to subset boundary geometries).
@@ -769,7 +773,6 @@
b = npy.asarray(self._boundarypolyll.boundary)
blons = b[:,0]; blats = b[:,1]
boundarypolyxy = PolygonShape(zip(*maptran(blons,blats)))
- # iterate over boundary geometries.
for line in bdatmetafile:
linesplit = line.split()
area = float(linesplit[1])
@@ -794,7 +797,7 @@
if not containsPole:
# close Antarctica for projections in which this
# is necessary.
- if name == 'gshhs' and self.projection in ['cyl','merc','mill','moll','robin','sinu','geos']:
+ if name == 'gshhs' and self.projection in ['cyl','merc','mill','moll','robin','sinu','geos','tmerc','cass','omerc']:
b = npy.asarray(poly.boundary)
lons = b[:,0]
lats = b[:,1]
@@ -817,8 +820,8 @@
continue
# if polygon instersects map projection
# region, process it.
- if poly.intersects(self._boundarypolyll):
- poly = poly.intersection(self._boundarypolyll)
+ if poly.intersects(boundarypolyll):
+ poly = poly.intersection(boundarypolyll)
# create iterable object with geometries
# that intersect map region.
if hasattr(poly,'geoms'):
@@ -887,17 +890,15 @@
continue
# create a Line object for other geometries.
poly = LineShape(zip(bx,by))
- # if polygon instersects map projection
- # region, process it.
- 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)
+ # if geometry instersects map projection
+ # region, and doesn't have any invalid points, process it.
+ if not badmask.any() and boundarypolyxy.intersects(poly):
+ # if geometry intersection calculation fails,
+ # just skip this polygon.
+ try:
+ poly = boundarypolyxy.intersection(poly)
+ except:
+ continue
# create iterable object with geometries
# that intersect map region.
if hasattr(poly,'geoms'):
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-07 13:11:21
|
Revision: 4140
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4140&view=rev
Author: jswhit
Date: 2007-11-07 05:11:15 -0800 (Wed, 07 Nov 2007)
Log Message:
-----------
add to list of projections that can't cross pole,
close Antarctica for all that don't.
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-07 12:56:23 UTC (rev 4139)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-07 13:11:15 UTC (rev 4140)
@@ -754,12 +754,9 @@
hasSP = boundarypolyxy.contains(SPole)
containsPole = hasNP or hasSP
# these projections cannot cross pole.
- if containsPole and self.projection in ['tmerc','cass','omerc','merc']:
+ if containsPole and\
+ self.projection in ['tmerc','cass','omerc','merc','mill','cyl','robin','moll','sinu','geos']:
raise ValueError('%s projection cannot cross pole'%(self.projection))
- # make sure geostationary projection has containsPole=False
- # even if only a limited area is requested (since entire
- # disk will be used to subset boundary geometries).
- if self.projection == 'geos': containsPole = False
# make sure orthographic projection has containsPole=True
# we will compute the intersections in stereographic
# coordinates, then transform to orthographic.
@@ -795,9 +792,9 @@
# coordinates (this saves time, especially for small map
# regions and high-resolution boundary geometries).
if not containsPole:
- # close Antarctica for projections in which this
+ # close Antarctica. for projections in which this
# is necessary.
- if name == 'gshhs' and self.projection in ['cyl','merc','mill','moll','robin','sinu','geos','tmerc','cass','omerc']:
+ if name == 'gshhs':
b = npy.asarray(poly.boundary)
lons = b[:,0]
lats = b[:,1]
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-07 16:13:00
|
Revision: 4143
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4143&view=rev
Author: jswhit
Date: 2007-11-07 08:12:57 -0800 (Wed, 07 Nov 2007)
Log Message:
-----------
optimize code for closing Antarctica
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-07 15:31:37 UTC (rev 4142)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-07 16:12:57 UTC (rev 4143)
@@ -18,6 +18,7 @@
from shapely.geometry import LineString as LineShape
from shapely.geometry import Point as PointShape
from shapely import wkb
+import time
# basemap data files now installed in lib/matplotlib/toolkits/basemap/data
basemap_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
@@ -743,6 +744,7 @@
bdatmetafile = open(os.path.join(basemap_datadir,name+'meta_'+self.resolution+'.dat'),'r')
except:
raise IOError, msg
+ timetot = 0.
polygons = []
polygon_types = []
# see if map projection region polygon contains a pole.
@@ -792,29 +794,25 @@
# coordinates (this saves time, especially for small map
# regions and high-resolution boundary geometries).
if not containsPole:
- # close Antarctica. for projections in which this
- # is necessary.
- if name == 'gshhs':
+ # close Antarctica.
+ if name == 'gshhs' and south < -68:
b = npy.asarray(poly.boundary)
lons = b[:,0]
lats = b[:,1]
- if south < -68:
- if math.fabs(lons[0]+0.) < 1.e-5:
- lons1 = lons[:-2][::-1]
- lats1 = lats[:-2][::-1]
- lons2 = lons1 + 360.
- lons3 = lons2 + 360.
- lons = lons1.tolist()+lons2.tolist()+lons3.tolist()
- lats = lats1.tolist()+lats1.tolist()+lats1.tolist()
- lonstart,latstart = lons[0], lats[0]
- lonend,latend = lons[-1], lats[-1]
- lons.insert(0,lonstart)
- lats.insert(0,-90.)
- lons.append(lonend)
- lats.append(-90.)
- poly = PolygonShape(zip(lons,lats))
- else:
- continue
+ if math.fabs(lons[0]+0.) < 1.e-5:
+ lons1 = lons[:-2][::-1]
+ lats1 = lats[:-2][::-1]
+ lons2 = lons1 + 360.
+ lons3 = lons2 + 360.
+ lons = lons1.tolist()+lons2.tolist()+lons3.tolist()
+ lats = lats1.tolist()+lats1.tolist()+lats1.tolist()
+ lonstart,latstart = lons[0], lats[0]
+ lonend,latend = lons[-1], lats[-1]
+ lons.insert(0,lonstart)
+ lats.insert(0,-90.)
+ lons.append(lonend)
+ lats.append(-90.)
+ poly = PolygonShape(zip(lons,lats))
# if polygon instersects map projection
# region, process it.
if poly.intersects(boundarypolyll):
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2007-11-07 16:26:10
|
Revision: 4144
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4144&view=rev
Author: jswhit
Date: 2007-11-07 08:25:58 -0800 (Wed, 07 Nov 2007)
Log Message:
-----------
fix Antarctic polygon closing code
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-07 16:12:57 UTC (rev 4143)
+++ trunk/toolkits/basemap-testing/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-07 16:25:58 UTC (rev 4144)
@@ -813,6 +813,8 @@
lons.append(lonend)
lats.append(-90.)
poly = PolygonShape(zip(lons,lats))
+ else:
+ continue
# if polygon instersects map projection
# region, process it.
if poly.intersects(boundarypolyll):
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X