plotsst.py 1.79 KB
Newer Older
1 2
from __future__ import (absolute_import, division, print_function)

3
from mpl_toolkits.basemap import Basemap
4
from netCDF4 import Dataset, date2index
5 6
import numpy as np
import matplotlib.pyplot as plt
7 8 9
from datetime import datetime
date = datetime(2007,12,15,0) # date to plot.
# open dataset.
10
dataset = \
11
Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg')
12 13
timevar = dataset.variables['time']
timeindex = date2index(date,timevar) # find time index for desired date.
14
# read sst.  Will automatically create a masked array using
15
# missing_value variable attribute. 'squeeze out' singleton dimensions.
16
sst = dataset.variables['sst'][timeindex,:].squeeze()
17
# read ice.
18
ice = dataset.variables['ice'][timeindex,:].squeeze()
19 20 21
# read lats and lons (representing centers of grid boxes).
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
22
lons, lats = np.meshgrid(lons,lats)
23
# create figure, axes instances.
24 25
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
26
# create Basemap instance.
27 28
# coastlines not used, so resolution set to None to skip
# continent processing (this speeds things up a bit)
29
m = Basemap(projection='kav7',lon_0=0,resolution=None)
30 31 32 33 34
# draw line around map projection limb.
# color background of map projection region.
# missing values over land will show up this color.
m.drawmapboundary(fill_color='0.3')
# plot sst, then ice with pcolor
35 36
im1 = m.pcolormesh(lons,lats,sst,shading='flat',cmap=plt.cm.jet,latlon=True)
im2 = m.pcolormesh(lons,lats,ice,shading='flat',cmap=plt.cm.gist_gray,latlon=True)
37
# draw parallels and meridians, but don't bother labelling them.
38 39
m.drawparallels(np.arange(-90.,99.,30.))
m.drawmeridians(np.arange(-180.,180.,60.))
40 41 42
# add colorbar
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
# add a title.
43
ax.set_title('SST and ICE analysis for %s'%date)
44
plt.show()