Menu
▾
▴
SF.net SVN: matplotlib:[8125] trunk/toolkits/basemap
|
From: <js...@us...> - 2010-02-11 13:03:03
|
Revision: 8125
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=8125&view=rev
Author: jswhit
Date: 2010-02-11 13:02:50 +0000 (Thu, 11 Feb 2010)
Log Message:
-----------
initial support for near-sided perspective satellite projection.
Modified Paths:
--------------
trunk/toolkits/basemap/Changelog
trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
trunk/toolkits/basemap/lib/mpl_toolkits/basemap/proj.py
Added Paths:
-----------
trunk/toolkits/basemap/examples/nsper_demo.py
Modified: trunk/toolkits/basemap/Changelog
===================================================================
--- trunk/toolkits/basemap/Changelog 2010-02-10 16:00:15 UTC (rev 8124)
+++ trunk/toolkits/basemap/Changelog 2010-02-11 13:02:50 UTC (rev 8125)
@@ -1,4 +1,6 @@
version 0.99.5 (not yet released)
+ * added "near-sided perspective" projection for a satellite
+ view at an arbitrary altitude.
* patch from Stephane Raynaud to pass format string to
drawmapscale, and allow units='m'.
* updated proj4 source to version 4.7.0, pyproj to 1.8.6.
Added: trunk/toolkits/basemap/examples/nsper_demo.py
===================================================================
--- trunk/toolkits/basemap/examples/nsper_demo.py (rev 0)
+++ trunk/toolkits/basemap/examples/nsper_demo.py 2010-02-11 13:02:50 UTC (rev 8125)
@@ -0,0 +1,25 @@
+from mpl_toolkits.basemap import Basemap
+import numpy as np
+import matplotlib.pyplot as plt
+
+# create Basemap instance for Near-Sided Perspective (satellite view) projection.
+lon_0 = float(raw_input('enter reference longitude (lon_0):'))
+lat_0 = float(raw_input('enter reference longitude (lat_0):'))
+h = float(raw_input('enter altitude of camera in km (h):'))
+h=h*1000.
+
+# map with continents drawn and filled.
+fig = plt.figure()
+m = Basemap(projection='nsper',lon_0=lon_0,lat_0=lat_0,satellite_height=h,resolution='l')
+m.drawcoastlines()
+m.drawmapboundary(fill_color='aqua')
+m.fillcontinents(color='coral',lake_color='aqua')
+m.drawcountries()
+m.drawstates()
+# draw parallels and meridians.
+m.drawparallels(np.arange(-90.,120.,10.))
+m.drawmeridians(np.arange(0.,420.,20.))
+m.drawmapboundary()
+plt.title('Near-Sided Perspective Map Centered on Lon=%s, Lat=%s, H=%g' %\
+ (lon_0,lat_0,h),fontsize=10)
+plt.show()
Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py 2010-02-10 16:00:15 UTC (rev 8124)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py 2010-02-11 13:02:50 UTC (rev 8125)
@@ -78,6 +78,7 @@
'poly' : 'Polyconic',
'ortho' : 'Orthographic',
'geos' : 'Geostationary',
+ 'nsper' : 'Near-Sided Perspective',
'sinu' : 'Sinusoidal',
'moll' : 'Mollweide',
'robin' : 'Robinson',
@@ -116,6 +117,7 @@
'poly' : 'lon_0,lat_0',
'ortho' : 'lon_0,lat_0,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
'geos' : 'lon_0,satellite_height,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
+ 'nsper' : 'lon_0,satellite_height,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
'sinu' : 'lon_0,lat_0,no corners or width/height',
'moll' : 'lon_0,lat_0,no corners or width/height',
'robin' : 'lon_0,lat_0,no corners or width/height',
@@ -189,10 +191,10 @@
For the cylindrical projections (``cyl``, ``merc``, ``mill`` and ``gall``),
the default is to use
llcrnrlon=-180,llcrnrlat=-90, urcrnrlon=180 and urcrnrlat=90). For all other
- projections except ``ortho`` and ``geos``, either the lat/lon values of the
+ projections except ``ortho``, ``geos`` and ``nsper``, either the lat/lon values of the
corners or width and height must be specified by the user.
- For ``ortho`` and ``geos``, the lat/lon values of the corners may be specified,
+ For ``ortho``, ``geos`` and ``nsper``, the lat/lon values of the corners may be specified,
or the x/y values of the corners (llcrnrx,llcrnry,urcrnrx,urcrnry) in the
coordinate system of the global projection (with x=0,y=0 at the center
of the global projection). If the corners are not specified,
@@ -311,8 +313,9 @@
latitude circle boundinglat is tangent to the edge
of the map at lon_0.
satellite_height height of satellite (in m) above equator -
- only relevant for geostationary projections
- (``geos``). Default 35,786 km.
+ only relevant for geostationary
+ and near-sided perspective (``geos`` or ``nsper``)
+ projections. Default 35,786 km.
================ ====================================================
Useful instance variables:
@@ -475,7 +478,7 @@
_insert_validated(projparams, lon_0, 'lon_0', -360, 720)
_insert_validated(projparams, lon_1, 'lon_1', -360, 720)
_insert_validated(projparams, lon_2, 'lon_2', -360, 720)
- if projection == 'geos':
+ if projection in ['geos','nsper']:
projparams['h'] = satellite_height
# check for sane values of projection corners.
using_corners = (None not in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat])
@@ -601,6 +604,24 @@
self._fulldisk = False
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
+ elif projection == 'nsper':
+ if not projparams.has_key('R'):
+ raise ValueError, 'near-sided perspective projection only works for perfect spheres - not ellipsoids'
+ if lat_0 is None or lon_0 is None:
+ msg='must specify lon_0 and lat_0 for near-sided perspective Basemap'
+ raise ValueError(msg)
+ if width is not None or height is not None:
+ print 'warning: width and height keywords ignored for %s projection' % _projnames[self.projection]
+ if not using_corners:
+ llcrnrlon = -180.
+ llcrnrlat = -90.
+ urcrnrlon = 180
+ urcrnrlat = 90.
+ self._fulldisk = True
+ else:
+ self._fulldisk = False
+ self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
+ self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection in _pseudocyl:
if lon_0 is None:
raise ValueError, 'must specify lon_0 for %s projection' % _projnames[self.projection]
@@ -723,7 +744,7 @@
self.aspect = (self.urcrnrlat-self.llcrnrlat)/(self.urcrnrlon-self.llcrnrlon)
else:
self.aspect = (proj.ymax-proj.ymin)/(proj.xmax-proj.xmin)
- if projection in ['geos','ortho'] and \
+ if projection in ['geos','ortho','nsper'] and \
None not in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
self.llcrnrx = llcrnrx+0.5*proj.xmax
self.llcrnry = llcrnry+0.5*proj.ymax
@@ -764,7 +785,7 @@
self.latmax = self.urcrnrlat
self.lonmin = self.llcrnrlon
self.lonmax = self.urcrnrlon
- elif self.projection in ['ortho','geos'] + _pseudocyl:
+ elif self.projection in ['ortho','geos','nsper'] + _pseudocyl:
self.latmin = -90.
self.latmax = 90.
self.lonmin = self.llcrnrlon
@@ -906,12 +927,12 @@
if containsPole and\
self.projection in _cylproj + _pseudocyl + ['geos']:
raise ValueError('%s projection cannot cross pole'%(self.projection))
- # make sure orthographic or gnomonic projection has containsPole=True
+ # make sure some projections have has containsPole=True
# we will compute the intersections in stereographic
- # coordinates, then transform to orthographic. This is
+ # coordinates, then transform back. This is
# because these projections are only defined on a hemisphere, and
# some boundary features (like Eurasia) would be undefined otherwise.
- if self.projection in ['ortho','gnom'] and name == 'gshhs':
+ if self.projection in ['ortho','gnom','nsper'] and name == 'gshhs':
containsPole = True
lon_0=self.projparams['lon_0']
lat_0=self.projparams['lat_0']
@@ -922,7 +943,7 @@
lat0 = 90.*(np.around(lat_0/90.))
if np.abs(int(lat0)) == 90: lon0=0.
maptran = pyproj.Proj(proj='stere',lon_0=lon0,lat_0=lat0,R=re)
- # boundary polygon for ortho/gnom projection
+ # boundary polygon for ortho/gnom/nsper projection
# in stereographic coorindates.
b = self._boundarypolyll.boundary
blons = b[:,0]; blats = b[:,1]
@@ -1030,9 +1051,10 @@
else:
# transform coordinates from lat/lon
# to map projection coordinates.
- # special case for ortho/gnom, compute coastline polygon
+ # special case for ortho/gnom/nsper, compute coastline polygon
# vertices in stereographic coords.
- if name == 'gshhs' and self.projection in ['ortho','gnom']:
+ if name == 'gshhs' and self.projection in\
+ ['ortho','gnom','nsper']:
b[:,0], b[:,1] = maptran(b[:,0], b[:,1])
else:
b[:,0], b[:,1] = self(b[:,0], b[:,1])
@@ -1046,10 +1068,10 @@
# in this map projection.
bx = np.compress(goodmask, b[:,0])
by = np.compress(goodmask, b[:,1])
- # for ortho/gnom projection, all points
+ # for ortho/gnom/nsper projection, all points
# outside map projection region are eliminated
# with the above step, so we're done.
- if self.projection in ['ortho','gnom']:
+ if self.projection in ['ortho','gnom','nsper']:
polygons.append(zip(bx,by))
polygon_types.append(type)
continue
@@ -1073,22 +1095,22 @@
# iterate over geometries in intersection.
for psub in geoms:
b = psub.boundary
- # if projection in ['ortho','gnom'],
+ # if projection in ['ortho','gnom','nsper'],
# transform polygon from stereographic
- # to ortho/gnom coordinates.
- if self.projection in ['ortho','gnom']:
+ # to ortho/gnom/nsper coordinates.
+ if self.projection in ['ortho','gnom','nsper']:
# if coastline polygon covers more than 99%
# of map region for fulldisk projection,
# it's probably bogus, so skip it.
areafrac = psub.area()/boundarypolyxy.area()
- if self.projection == 'ortho':
+ if self.projection == ['ortho','nsper']:
if name == 'gshhs' and\
self._fulldisk and\
areafrac > 0.99: continue
# inverse transform from stereographic
# to lat/lon.
b[:,0], b[:,1] = maptran(b[:,0], b[:,1], inverse=True)
- # orthographic/gnomonic.
+ # orthographic/gnomonic/nsper.
b[:,0], b[:,1]= self(b[:,0], b[:,1])
polygons.append(zip(b[:,0],b[:,1]))
polygon_types.append(type)
@@ -1102,7 +1124,7 @@
if self.projection == 'vandg':
nx = 10*nx; ny = 10*ny
maptran = self
- if self.projection in ['ortho','geos']:
+ if self.projection in ['ortho','geos','nsper']:
# circular region.
thetas = np.linspace(0.,2.*np.pi,2*nx*ny)[:-1]
rminor = self._height
@@ -1246,10 +1268,10 @@
# get current axes instance (if none specified).
ax = ax or self._check_ax()
limb = None
- if self.projection in ['ortho','geos'] or (self.projection=='aeqd' and\
+ if self.projection in ['ortho','geos','nsper'] or (self.projection=='aeqd' and\
self._fulldisk):
limb = Ellipse((self._width,self._height),2.*self._width,2.*self._height)
- if self.projection in ['ortho','geos','aeqd'] and self._fulldisk:
+ if self.projection in ['ortho','geos','nsper','aeqd'] and self._fulldisk:
# elliptical region.
ax.add_patch(limb)
if fill_color is None:
@@ -1301,7 +1323,7 @@
except AttributeError:
for spine in ax.spines.itervalues():
spine.set_linewidth(linewidth)
- if self.projection not in ['geos','ortho']:
+ if self.projection not in ['geos','ortho','nsper']:
if fill_color is not None:
ax.axesPatch.set_facecolor(fill_color)
try:
@@ -1876,7 +1898,7 @@
linecolls[circ] = (lines,[])
# draw labels for parallels
# parallels not labelled for fulldisk orthographic or geostationary
- if self.projection in ['ortho','geos','vandg','aeqd'] and max(labels):
+ if self.projection in ['ortho','geos','nsper','vandg','aeqd'] and max(labels):
if self.projection == 'vandg' or self._fulldisk:
print 'Warning: Cannot label parallels on %s basemap' % _projnames[self.projection]
labels = [0,0,0,0]
@@ -2118,7 +2140,7 @@
if self.projection in ['sinu','moll','vandg'] and max(labels):
print 'Warning: Cannot label meridians on %s basemap' % _projnames[self.projection]
labels = [0,0,0,0]
- if self.projection in ['ortho','geos','aeqd'] and max(labels):
+ if self.projection in ['ortho','geos','nsper','aeqd'] and max(labels):
if self._fulldisk:
print dedent(
"""'Warning: Cannot label meridians on full-disk
@@ -2572,7 +2594,7 @@
# turn off axes frame for non-rectangular projections.
if self.projection in _pseudocyl:
ax.set_frame_on(False)
- if self.projection in ['ortho','geos','aeqd'] and self._fulldisk:
+ if self.projection in ['ortho','geos','nsper','aeqd'] and self._fulldisk:
ax.set_frame_on(False)
# make sure aspect ratio of map preserved.
# plot is re-centered in bounding rectangle.
@@ -3298,7 +3320,7 @@
self.transform_scalar(self._bm_rgba[:,:,k],\
self._bm_lons,self._bm_lats,nx,ny,returnxy=True)
# for ortho,geos mask pixels outside projection limb.
- if self.projection in ['geos','ortho'] or \
+ if self.projection in ['geos','ortho','nsper'] or \
(self.projection == 'aeqd' and self._fulldisk):
lonsr,latsr = self(x,y,inverse=True)
mask = ma.zeros((ny,nx,4),np.int8)
Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/proj.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/proj.py 2010-02-10 16:00:15 UTC (rev 8124)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/proj.py 2010-02-11 13:02:50 UTC (rev 8125)
@@ -147,6 +147,38 @@
llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
if llcrnrx > 1.e20 or llcrnry > 1.e20:
raise ValueError(_lower_left_out_of_bounds)
+ elif self.projection == 'nsper':
+ self._proj4 = pyproj.Proj(projparams)
+ # find major and minor axes of ellipse defining map proj region.
+ # h is measured from surface of earth at equator.
+ h = projparams['h'] + self.rmajor
+ # latitude of horizon on central meridian
+ lonmax = 90.-(180./np.pi)*np.arcsin(self.rmajor/h)
+ # longitude of horizon on equator
+ latmax = 90.-(180./np.pi)*np.arcsin(self.rmajor/h)
+ # truncate to nearest hundredth of a degree (to make sure
+ # they aren't slightly over the horizon)
+ latmax = int(100*latmax)/100.
+ lonmax = int(100*lonmax)/100.
+ # width and height of visible projection
+ P = pyproj.Proj(proj='nsper',a=self.rmajor,\
+ b=self.rminor,lat_0=0,lon_0=0,h=projparams['h'])
+ x1,y1 = P(0.,latmax); x2,y2 = P(lonmax,0.)
+ width = x2; height = y1
+ self._height = height
+ self._width = width
+ if (llcrnrlon == -180 and llcrnrlat == -90 and
+ urcrnrlon == 180 and urcrnrlat == 90):
+ self._fulldisk = True
+ llcrnrx = -width
+ llcrnry = -height
+ urcrnrx = -llcrnrx
+ urcrnry = -llcrnry
+ else:
+ self._fulldisk = False
+ llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
+ if llcrnrx > 1.e20 or llcrnry > 1.e20:
+ raise ValueError(_lower_left_out_of_bounds)
elif self.projection in _pseudocyl:
self._proj4 = pyproj.Proj(projparams)
xtmp,urcrnry = self(projparams['lon_0'],90.)
@@ -172,9 +204,9 @@
if urcrnrislatlon:
self.urcrnrlon = urcrnrlon
self.urcrnrlat = urcrnrlat
- if self.projection not in ['ortho','geos','aeqd'] + _pseudocyl:
+ if self.projection not in ['ortho','geos','nsper','aeqd'] + _pseudocyl:
urcrnrx,urcrnry = self(urcrnrlon,urcrnrlat)
- elif self.projection in ['ortho','geos','aeqd']:
+ elif self.projection in ['ortho','geos','nsper','aeqd']:
if self._fulldisk:
urcrnrx = 2.*self._width
urcrnry = 2.*self._height
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
