Merge pull request #23 from jswhit/master

update coastlines to gshhs 2.2.0

Author: Jeff Whitaker

Changelog

@@ -1,13 +1,17 @@
 version 1.0.3 (not yet released)
-           * added hexbin method.
-           * drawmapboundary now uses axes bgcolor as default fill_color. If
-             no color fill wanted, set fill_color='none' (a string).
-           * clip coastlines for nplaea,npaeqd,splaea,spaeqd in stereographic
-             coordinates to avoid S. America disappearing in some south polar
-             plots.
-           * add 'round' keyword to Basemap.__init__ for pole-centered
+           * update coastlines, rivers, political boundaries to GSHHS
+             2.2.0/GMT 4.5.7. The fillcontinents bug (filling the outside
+	     instead of the inside of a coastline polygon) is now much
+             harder to trigger.    
+           * add 'round' keyword keyword to Basemap.__init__ for pole-centered
              projections to make them round (clipped at boundinglat) instead
-             of square.
+             of square.        
+	   * added hexbin method.
+	   * drawmapboundary now uses axes bgcolor as default fill_color. If
+	     no color fill wanted, set fill_color='none' (a string).
+	   * clip coastlines for nplaea,npaeqd,splaea,spaeqd in stereographic
+	     coordinates to avoid S. America disappearing in some south polar
+	     plots.
            * fix broken daynight terminator function.
            * added kav7 (Kavrayskiy VII) and eck4 (Eckert IV) projections
              (https://github.com/matplotlib/basemap/issues/9).

KNOWN_BUGS

@@ -1,7 +1,6 @@
 * The fillcontinents method doesn't always do the right thing.  Matplotlib
 always tries to fill the inside of a polygon.  Under certain situations,
 what is the inside of a coastline polygon can be ambiguous, and the 
-outside may be filled instead of the inside.  To trigger this,
-run the ortho_demo.py example with lon=-120,lat=60.
+outside may be filled instead of the inside. 
 Workaround - mask the land areas with the drawlsmask method instead of
 filling the coastline polygons (this is illustrated in the ortho_demo.py example). 

doc/users/figures/plotsst.py

@@ -1,5 +1,5 @@
 from mpl_toolkits.basemap import Basemap
-from netCDF4 import Dataset, date2index, num2date
+from netCDF4 import Dataset
 import numpy as np
 import matplotlib.pyplot as plt
 date = '20071215' # date to plot.

doc/users/figures/plotwindvec.py

@@ -2,33 +2,31 @@
 import matplotlib.pyplot as plt
 import datetime
 from mpl_toolkits.basemap import Basemap, shiftgrid
-from netCDF4 import Dataset, date2index, num2date
+from netCDF4 import Dataset
 # specify date to plot.
-date = datetime.datetime(1993,3,14,0)
+yyyy=1993; mm=03; dd=14; hh=00
+date = datetime.datetime(yyyy,mm,dd,hh)
 # set OpenDAP server URL.
-URL="http://nomad1.ncep.noaa.gov:9090/dods/reanalyses/reanalysis-2/6hr/pgb/pgb"
+URLbase="http://nomads.ncdc.noaa.gov/thredds/dodsC/modeldata/cmd_pgbh/"
+URL=URLbase+"%04i/%04i%02i/%04i%02i%02i/pgbh00.gdas.%04i%02i%02i%02i.grb2" %\
+             (yyyy,yyyy,mm,yyyy,mm,dd,yyyy,mm,dd,hh)
 data = Dataset(URL)
-# read lats,lons,times.
-latitudes = data.variables['lat'][:]
+# read lats,lons
+# reverse latitudes so they go from south to north.
+latitudes = data.variables['lat'][::-1]
 longitudes = data.variables['lon'][:].tolist()
-times = data.variables['time']
-# get time index for desired date.
-ntime = date2index(date,times,calendar='standard')
 # get sea level pressure and 10-m wind data.
-slpdata = data.variables['presmsl']
-udata = data.variables['ugrdprs']
-vdata = data.variables['vgrdprs']
 # mult slp by 0.01 to put in units of hPa.
-slpin = 0.01*slpdata[ntime,:,:]
-uin = udata[ntime,0,:,:] 
-vin = vdata[ntime,0,:,:] 
+slpin = 0.01*data.variables['Pressure_msl'][:].squeeze()
+uin = data.variables['U-component_of_wind_height_above_ground'][:].squeeze()
+vin = data.variables['V-component_of_wind_height_above_ground'][:].squeeze()
 # add cyclic points manually (could use addcyclic function)
 slp = np.zeros((slpin.shape[0],slpin.shape[1]+1),np.float)
-slp[:,0:-1] = slpin; slp[:,-1] = slpin[:,0]
+slp[:,0:-1] = slpin[::-1]; slp[:,-1] = slpin[::-1,0]
 u = np.zeros((uin.shape[0],uin.shape[1]+1),np.float64)
-u[:,0:-1] = uin; u[:,-1] = uin[:,0]
+u[:,0:-1] = uin[::-1]; u[:,-1] = uin[::-1,0]
 v = np.zeros((vin.shape[0],vin.shape[1]+1),np.float64)
-v[:,0:-1] = vin; v[:,-1] = vin[:,0]
+v[:,0:-1] = vin[::-1]; v[:,-1] = vin[::-1,0]
 longitudes.append(360.); longitudes = np.array(longitudes)
 # make 2-d grid of lons, lats
 lons, lats = np.meshgrid(longitudes,latitudes)

examples/fcstmaps.py

@@ -59,13 +59,12 @@
 # unpack 2-meter temp forecast data.
 
 t2mvar = data.variables['tmp2m']
-t2min = t2mvar[0:ntimes,:,:]
-t2m = np.zeros((ntimes,nlats,nlons+1),t2min.dtype)
+t2m = np.zeros((ntimes,nlats,nlons+1),np.float32)
 # create Basemap instance for Orthographic projection.
 m = Basemap(lon_0=-90,lat_0=60,projection='ortho')
 # add wrap-around point in longitude.
 for nt in range(ntimes):
-    t2m[nt,:,:], lons = addcyclic(t2min[nt,:,:], lons1)
+    t2m[nt,:,:], lons = addcyclic(t2mvar[nt,:,:], lons1)
 # convert to celsius.
 t2m = t2m-273.15
 # contour levels

examples/fcstmaps_axesgrid.py

@@ -61,8 +61,6 @@
 # unpack 2-meter temp forecast data.
 
 t2mvar = data.variables['tmp2m']
-t2min = t2mvar[0:ntimes,:,:]
-t2m = np.zeros((ntimes,nlats,nlons+1),t2min.dtype)
 
 # create figure, set up AxesGrid.
 fig=plt.figure(figsize=(6,8))
@@ -79,8 +77,9 @@
 # create Basemap instance for Orthographic projection.
 m = Basemap(lon_0=-90,lat_0=60,projection='ortho')
 # add wrap-around point in longitude.
+t2m = np.zeros((ntimes,nlats,nlons+1),np.float32)
 for nt in range(ntimes):
-    t2m[nt,:,:], lons = addcyclic(t2min[nt,:,:], lons1)
+    t2m[nt,:,:], lons = addcyclic(t2mvar[nt,:,:], lons1)
 # convert to celsius.
 t2m = t2m-273.15
 # contour levels

examples/nsper_demo.py

@@ -4,7 +4,7 @@
 
 # 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):'))
+lat_0 = float(raw_input('enter reference latitude (lat_0):'))
 h = float(raw_input('enter altitude of camera in km (h):'))
 h=h*1000.
 

examples/plotsst.py

@@ -1,21 +1,17 @@
 from mpl_toolkits.basemap import Basemap
-from netCDF4 import Dataset, date2index, num2date
+from netCDF4 import Dataset
 import numpy as np
 import matplotlib.pyplot as plt
-import datetime
-# create datetime object for desired time
-date = datetime.datetime(2007,12,15,0)
-# open dataset.
-dataset =\
-Dataset('http://nomads.ncdc.noaa.gov/thredds/dodsC/oisst2/totalAmsrAgg')
-# find index of desired time.
-time = dataset.variables['time']
-nt = date2index(date, time, calendar='standard')
+date = '20071215' # date to plot.
+# open dataset for that date.
+dataset = \
+Dataset('http://nomads.ncdc.noaa.gov/thredds/dodsC/oisst2/%s/AVHRR-AMSR/amsr-avhrr-v2.%s.nc'%\
+       (date[0:4],date))
 # read sst.  Will automatically create a masked array using
-# missing_value variable attribute.
-sst = dataset.variables['sst'][nt]
+# missing_value variable attribute. 'squeeze out' singleton dimensions.
+sst = dataset.variables['sst'][:].squeeze()
 # read ice.
-ice = dataset.variables['ice'][nt]
+ice = dataset.variables['ice'][:].squeeze()
 # read lats and lons (representing centers of grid boxes).
 lats = dataset.variables['lat'][:]
 lons = dataset.variables['lon'][:]
@@ -28,7 +24,7 @@
 lats = (lats - 0.5*delat).tolist()
 lats.append(lats[-1]+delat)
 lats = np.array(lats,np.float64)
-# creat figure, axes instances.
+# create figure, axes instances.
 fig = plt.figure()
 ax = fig.add_axes([0.05,0.05,0.9,0.9])
 # create Basemap instance for Robinson projection.
@@ -50,5 +46,5 @@
 # add colorbar
 cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
 # add a title.
-plt.title('SST and ICE analysis for %s'%date)
+ax.set_title('SST and ICE analysis for %s'%date)
 plt.show()