Merge pull request #105 from jswhit/master

update test_rotpole.py example

Author: Jeff Whitaker

examples/test_rotpole.py

@@ -1,3 +1,4 @@
+from __future__ import print_function
 from netCDF4 import Dataset
 from mpl_toolkits.basemap import Basemap
 import numpy as np
@@ -23,17 +24,33 @@
 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].
+#  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')
+
+def normalize180(lon):
+    """Normalize lon to range [180, 180)"""
+    lower = -180.; upper = 180.
+    if lon > upper or lon == lower:
+        lon = lower + abs(lon + upper) % (abs(lower) + abs(upper))
+    if lon < lower or lon == upper:
+        lon = upper - abs(lon - lower) % (abs(lower) + abs(upper))
+    return lower if lon == upper else lon
+
+lon_0 = normalize180(rotpole.grid_north_pole_longitude-180.)
+o_lon_p = rotpole.north_pole_grid_longitude
+o_lat_p = rotpole.grid_north_pole_latitude
+print( rotpole )
+print( 'lon_0,o_lon_p,o_lat_p=',lon_0,o_lon_p,o_lat_p)
+
+m= Basemap(projection='rotpole',lon_0=lon_0,o_lon_p=o_lon_p,o_lat_p=o_lat_p,\
+           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)
@@ -43,13 +60,9 @@
 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
+m= Basemap(projection='rotpole',lon_0=lon_0,o_lon_p=o_lon_p,o_lat_p=o_lat_p,\
+           llcrnry = rlats[0,0], urcrnry = rlats[-1,-1],\
+           llcrnrx = rlons[0,0], urcrnrx = rlons[-1,-1],resolution='c')
 x,y = m(lons,lats)
 m.drawcoastlines()
 m.contourf(x,y,data,20)