Merge pull request #104 from jswhit/master

add projection='rotpole' (rotated pole transformation)

Author: Jeff Whitaker

Changelog

@@ -1,5 +1,6 @@
 since version 1.0.6
 -------------------
+* support for rotated pole transormation (projection = 'rotpole').
 * fix warpimage with projection = 'hammer' (issue 100).
 
 version 1.0.6 (git tag v1.0.6rel)

MANIFEST.in

@@ -88,6 +88,8 @@ include examples/ccsm_popgrid.nc
 include examples/rita.nc
 include examples/maskoceans.py
 include examples/utmtest.py
+include examples/test_rotpole.py
+include examples/wm201_Arctic_JJA_1990-2008_moyenneDesMoyennes.nc
 include examples/README
 include lib/mpl_toolkits/basemap/__init__.py
 include lib/mpl_toolkits/basemap/proj.py

examples/README

@@ -164,3 +164,5 @@ ArcGIS server using the REST API.
 
 testwmsimage.py shows how to retrieve an image from a WMS server and display
 it on a map.
+
+test_rotpole.py shows how to plot regional climate model data in the native 'rotated pole' projection.

examples/test.py

@@ -369,6 +369,22 @@
 print('plotting Square Lambert Azimuthal example ...')
 print(m.proj4string)
 
+# create new figure
+fig=plt.figure()
+m = Basemap(projection = 'rotpole',lon_0 = -120.,\
+            o_lon_p = 180, o_lat_p = 0,\
+            llcrnry = -41.75, urcrnry = 37.75,\
+            llcrnrx = 137, urcrnrx = 222.5, resolution = 'l')
+m.drawcoastlines()
+ny,nx = lons.shape
+m.contourf(lons[ny/2:,:],lats[ny/2:,:],topodat[ny/2:,:],50,cmap=cmap,extend='both',latlon=True)
+m.drawmeridians(np.arange(-180,180,20),labels=[1,1,1,1])
+m.drawparallels(np.arange(20,80,20))
+m.colorbar()
+plt.title('Rotated Pole',y=1.075)
+print('plotting Rotated Pole example ...')
+print(m.proj4string)
+
 # create new figure
 fig=plt.figure()
 m = Basemap(lon_0=-105,boundinglat=20.,

examples/test_rotpole.py

@@ -0,0 +1,61 @@
+from netCDF4 import Dataset
+from mpl_toolkits.basemap import Basemap
+import numpy as np
+import matplotlib.pyplot as plt
+
+nc = Dataset('wm201_Arctic_JJA_1990-2008_moyenneDesMoyennes.nc')
+lats = nc.variables['lat'][:]
+lons = nc.variables['lon'][:]
+rlats = nc.variables['rlat'][:]
+rlons = nc.variables['rlon'][:]
+rlons, rlats = np.meshgrid(rlons, rlats)
+data = nc.variables['air'][0,0,:,:].squeeze()
+data = np.ma.masked_values(data,-999.)
+rotpole = nc.variables['rotated_pole']
+
+m = Basemap(projection='npstere',lon_0=10,boundinglat=30,resolution='c')
+x,y = m(lons,lats)
+m.drawcoastlines()
+m.contourf(x,y,data,20)
+m.drawmeridians(np.arange(-180,180,20))
+m.drawparallels(np.arange(20,80,20))
+m.colorbar()
+plt.title('rotated pole data in polar stere map')
+
+plt.figure()
+# o_lon_p, o_lat_p: true lat/lon of pole in rotated coordinate system
+# mapping to CF metadata convention:
+# grid_north_pole_longitude = normalize180(180 + lon_0), where normalize180
+#  is a function that maps to interval [-180,180].
+# grid_north_pole_latitude = o_lat_p
+# north_pole_grid_longitude = o_lon_p (optional, assumed zero if not present)
+m = Basemap(projection='rotpole',lon_0=rotpole.grid_north_pole_longitude-180.,\
+            o_lon_p=rotpole.north_pole_grid_longitude,\
+            o_lat_p=rotpole.grid_north_pole_latitude,\
+            llcrnrlat = lats[0,0], urcrnrlat = lats[-1,-1],\
+            llcrnrlon = lons[0,0], urcrnrlon = lons[-1,-1],resolution='c')
+x,y = m(lons,lats)
+m.drawcoastlines()
+m.contourf(x,y,data,20)
+m.drawmeridians(np.arange(-180,180,20))
+m.drawparallels(np.arange(20,80,20))
+m.colorbar()
+plt.title('rotated pole data in native map using real sphere corner lat/lons' )
+
+plt.figure()
+m = Basemap(projection='rotpole',lon_0=rotpole.grid_north_pole_longitude-180.,\
+            o_lon_p=rotpole.north_pole_grid_longitude,\
+            o_lat_p=rotpole.grid_north_pole_latitude,\
+            llcrnry = rlats[0,0], urcrnry = rlats[-1,-1],\
+            llcrnrx = rlons[0,0], urcrnrx = rlons[-1,-1],resolution='c')
+print m.llcrnrx,m.llcrnry
+print m.urcrnrx,m.urcrnry
+x,y = m(lons,lats)
+m.drawcoastlines()
+m.contourf(x,y,data,20)
+m.drawmeridians(np.arange(-180,180,20))
+m.drawparallels(np.arange(20,80,20))
+m.colorbar()
+plt.title('rotated pole data in native map using rotated sphere corner lat/lons' )
+
+plt.show()

examples/wm201_Arctic_JJA_1990-2008_moyenneDesMoyennes.nc

@@ -88,6 +88,7 @@
              'vandg'    : 'van der Grinten',
              'mbtfpq'   : 'McBryde-Thomas Flat-Polar Quartic',
              'gnom'     : 'Gnomonic',
+             'rotpole'  : 'Rotated Pole',
              }
 supported_projections = []
 for _items in _projnames.items():
@@ -96,6 +97,8 @@
 
 _cylproj = ['cyl','merc','mill','gall','cea']
 _pseudocyl = ['moll','robin','eck4','kav7','sinu','mbtfpq','vandg','hammer']
+_dg2rad = math.radians(1.)
+_rad2dg = math.degrees(1.)
 
 # projection specific parameters.
 projection_params = {'cyl'      : 'corners only (no width/height)',
@@ -131,6 +134,7 @@
              'vandg'    : 'lon_0,lat_0,no corners or width/height',
              'mbtfpq'   : 'lon_0,lat_0,no corners or width/height',
              'gnom'     : 'lon_0,lat_0',
+             'rotpole'  : 'lon_0,o_lat_p,o_lon_p,corner lat/lon or corner x,y (no width/height)'
              }
 
 # create dictionary that maps epsg codes to Basemap kwargs.
@@ -258,6 +262,11 @@
  of the global projection).  If the corners are not specified,
  the entire globe is plotted.
 
+ For ``rotpole``, the lat/lon values of the corners on the unrotated sphere
+ may be provided as llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat, or the lat/lon
+ values of the corners on the rotated sphere can be given as
+ llcrnrx,llcrnry,urcrnrx,urcrnry.
+
  Other keyword arguments:
 
  .. tabularcolumns:: |l|L|
@@ -564,6 +573,7 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None,
                        lat_1=None, lat_2=None,
                        lat_0=None, lon_0=None,
                        lon_1=None, lon_2=None,
+                       o_lon_p=None, o_lat_p=None,
                        k_0=None,
                        no_rot=False,
                        suppress_ticks=True,
@@ -906,6 +916,27 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None,
                 projparams['lon_0'] = lon_0
             else:
                 projparams['lon_0']=0.5*(llcrnrlon+urcrnrlon)
+        elif projection == 'rotpole':
+            if lon_0 is None or o_lon_p is None or o_lat_p is None:
+                msg='must specify lon_0,o_lat_p,o_lon_p for rotated pole Basemap'
+                raise ValueError(msg)
+            if width is not None or height is not None:
+                sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
+            projparams['lon_0']=lon_0
+            projparams['o_lon_p']=o_lon_p
+            projparams['o_lat_p']=o_lat_p
+            projparams['o_proj']='longlat'
+            projparams['proj']='ob_tran'
+            if not using_corners and None in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
+                raise ValueError('must specify lat/lon values of corners in degrees')
+            if None not in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
+                p = pyproj.Proj(projparams)
+                llcrnrx = _dg2rad*llcrnrx; llcrnry = _dg2rad*llcrnry
+                urcrnrx = _dg2rad*urcrnrx; urcrnry = _dg2rad*urcrnry
+                llcrnrlon, llcrnrlat = p(llcrnrx,llcrnry,inverse=True)
+                urcrnrlon, urcrnrlat = p(urcrnrx,urcrnry,inverse=True)
+                self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
+                self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
         else:
             raise ValueError(_unsupported_projection % projection)
 
@@ -968,6 +999,12 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None,
             self.urcrnrx = proj.urcrnrx
             self.urcrnry = proj.urcrnry
 
+        if self.projection == 'rotpole':
+            lon0,lat0 = self(0.5*(self.llcrnrx + self.urcrnrx),\
+                             0.5*(self.llcrnry + self.urcrnry),\
+                             inverse=True)
+            self.projparams['lat_0']=lat0
+
         # if ax == None, pyplot.gca may be used.
         self.ax = ax
         self.lsmask = None
@@ -1097,12 +1134,32 @@ def __call__(self,x,y,inverse=False):
                 x = 2.*lon_0-x
             except TypeError:
                 x = [2*lon_0-xx for xx in x]
+        if self.projection == 'rotpole' and inverse:
+            try:
+                x = _dg2rad*x
+            except TypeError:
+                x = [_dg2rad*xx for xx in x]
+            try:
+                y = _dg2rad*y
+            except TypeError:
+                y = [_dg2rad*yy for yy in y]
         xout,yout = self.projtran(x,y,inverse=inverse)
         if self.celestial and inverse:
             try:
                 xout = -2.*lon_0-xout
             except:
                 xout = [-2.*lon_0-xx for xx in xout]
+        if self.projection == 'rotpole' and not inverse:
+            try:
+                xout = _rad2dg*xout
+                xout = np.where(xout < 0., xout+360, xout)
+            except TypeError:
+                xout = [_rad2dg*xx for xx in xout]
+                xout = [xx+360. if xx < 0 else xx for xx in xout]
+            try:
+                yout = _rad2dg*yout
+            except TypeError:
+                yout = [_rad2dg*yy for yy in yout]
         return xout,yout
 
     def makegrid(self,nx,ny,returnxy=False):
@@ -1155,7 +1212,8 @@ def _readboundarydata(self,name,as_polygons=False):
         # 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.
-        tostere = ['omerc','ortho','gnom','nsper','nplaea','npaeqd','splaea','spaeqd']
+        tostere =\
+        ['omerc','ortho','gnom','nsper','nplaea','npaeqd','splaea','spaeqd']
         if self.projection in tostere and name == 'gshhs':
             containsPole = True
             lon_0=self.projparams['lon_0']
@@ -1336,7 +1394,7 @@ def _readboundarydata(self,name,as_polygons=False):
                         xd = (bx[1:]-bx[0:-1])**2
                         yd = (by[1:]-by[0:-1])**2
                         dist = np.sqrt(xd+yd)
-                        split = dist > 5000000.
+                        split = dist > 0.5*(self.xmax-self.xmin)
                         if np.sum(split) and self.projection not in _cylproj:
                             ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist()
                             iprev = 0

lib/mpl_toolkits/basemap/__init__.py

@@ -74,6 +74,11 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,
         if self.projection == 'cyl':
             llcrnrx = llcrnrlon
             llcrnry = llcrnrlat
+        elif self.projection == 'ob_tran':
+            self._proj4 = pyproj.Proj(projparams)
+            llcrnrx,llcrnry = self(llcrnrlon,llcrnrlat)
+            llcrnrx = _rad2dg*llcrnrx; llcrnry = _rad2dg*llcrnry
+            if llcrnrx < 0: llcrnrx = llcrnrx + 360
         elif self.projection in 'ortho':
             if (llcrnrlon == -180 and llcrnrlat == -90 and
                 urcrnrlon == 180 and urcrnrlat == 90):
@@ -194,21 +199,25 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,
             if self.projection == 'aeqd': self._fulldisk=False
         # compute x_0, y_0 so ll corner of domain is x=0,y=0.
         # note that for 'cyl' x,y == lon,lat
-        self.projparams['x_0']=-llcrnrx
-        self.projparams['y_0']=-llcrnry
+        if self.projection != 'ob_tran':
+            self.projparams['x_0']=-llcrnrx
+            self.projparams['y_0']=-llcrnry
         # reset with x_0, y_0.
-        if self.projection != 'cyl':
+        if self.projection not in ['cyl','ob_tran']:
             self._proj4 = pyproj.Proj(projparams)
             llcrnry = 0.
             llcrnrx = 0.
-        else:
+        elif self.projection != 'ob_tran':
             llcrnrx = llcrnrlon
             llcrnry = llcrnrlat
         if urcrnrislatlon:
             self.urcrnrlon = urcrnrlon
             self.urcrnrlat = urcrnrlat
             if self.projection not in ['ortho','geos','nsper','aeqd'] + _pseudocyl:
                 urcrnrx,urcrnry = self(urcrnrlon,urcrnrlat)
+                if self.projection == 'ob_tran':
+                    urcrnrx = _rad2dg*urcrnrx; urcrnry = _rad2dg*urcrnry
+                    if urcrnrx < 0: urcrnrx = urcrnrx + 360
             elif self.projection in ['ortho','geos','nsper','aeqd']:
                 if self._fulldisk:
                     urcrnrx = 2.*self._width
@@ -226,6 +235,7 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,
             urcrnrlon, urcrnrlat = self(urcrnrx, urcrnry, inverse=True)
             self.urcrnrlon = urcrnrlon
             self.urcrnrlat = urcrnrlat
+
         # corners of domain.
         self.llcrnrx = llcrnrx
         self.llcrnry = llcrnry