Examples of accessing Netcdf data via TRHEDDS/OPENDAP services in Python, and plotting in Basemaps
First, import libraries
Important note It looks like for users on Windows with Python 3.x this demo will not work. It will work on Windows with Python 2.7 however. If you are on Linux or Mac (or running 2.7 in Windows) you can add the Basemap package with conda using the command conda install Basemap.
In [ ]:
from mpl_toolkits.basemap import Basemap, shiftgrid
from netCDF4 import Dataset, date2index
import time
import sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from datetime import datetime, timedelta
%matplotlib notebook
For example, ocean data from NOAA' catalouge
In [ ]:
nc = Dataset('http://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst2Agg')
In [ ]:
nc.variables # See the metadata via this
In [ ]:
# Grab a slice of the SST to play with
# nc.variables['sst']
sst = nc.variables['sst'][0,:].squeeze()
# preview the data with an array plotting function from matplotlib
fig, ax = plt.subplots(1,1)
ax.imshow(np.flipud(sst))
In [ ]:
lon, lat = sst.shape
print("Number of (floating point value) pixels of AVHRR data retrieved: {0:10,}".format(lon * lat))
print("Size in memory: {0:3.1f} mb".format(16 * (lon*lat)/1000000)) # 16 bytes in a float, 1 million bytes in a megabyte
In [ ]:
%%time
# based on example at http://matplotlib.org/basemap/users/examples.html
date = datetime(2007,12,15,0) # Specify date to plot.
dataset = Dataset('http://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst2Agg')
timevar = dataset.variables['time']
timeindex = date2index(date, timevar) # find time index for desired date.
# read sst. Will automatically create a masked array using
# missing_value variable attribute. 'squeeze out' singleton dimensions.
sst = dataset.variables['sst'][timeindex,:].squeeze()
# read ice.
ice = dataset.variables['ice'][timeindex,:].squeeze()
# read lats and lons (representing centers of grid boxes).
lats = dataset.variables['latitude'][:]
lons = dataset.variables['longitude'][:]
lons, lats = np.meshgrid(lons,lats)
# create figure, axes instances.
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
# create Basemap instance.
# coastlines not used, so resolution set to None to skip
# continent processing (this speeds things up a bit)
m = Basemap(projection='kav7',lon_0=0,resolution=None)
m.drawmapboundary(fill_color='0.3') # color map background
# plot sst, then ice with pcolor
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)
# draw parallels and meridians, but don't bother labelling them.
#m.drawparallels(np.arange(-90.,99.,30.))
#m.drawmeridians(np.arange(-180.,180.,60.))
# add colorbar
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
# add a title.
ax.set_title('SST and ICE location on {0}'.format(date.date()))
plt.show()
In [ ]:
# specify date to plot.
date = datetime(1993, 3, 14, 0)
yyyy = date.year
mm = date.month
dd = date.day
hh = date.hour
# set OpenDAP server URL.
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)
latitudes = data.variables['lat'][::-1]
longitudes = data.variables['lon'][:].tolist()
# Get pressure and 10-m wind data
slpin = 0.01*data.variables['Pressure_msl'][:].squeeze() # 0.01* to convert to hPa
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[::-1]; slp[:,-1] = slpin[::-1,0]
u = np.zeros((uin.shape[0],uin.shape[1]+1),np.float64)
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[::-1]; v[:,-1] = vin[::-1,0]
longitudes.append(360.)
longitudes = np.array(longitudes)
lons, lats = np.meshgrid(longitudes,latitudes) # make 2-d grid of lons, lats
m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.)
# create figure, add axes
fig1 = plt.figure(figsize=(8,10))
ax = fig1.add_axes([0.1,0.1,0.8,0.8])
clevs = np.arange(960,1061,5)
x, y = m(lons, lats)
parallels = np.arange(-80.,90,20.)
meridians = np.arange(0.,360.,20.)
# plot SLP contours.
CS1 = m.contour(x,y,slp,clevs,linewidths=0.5,colors='k',animated=True)
CS2 = m.contourf(x,y,slp,clevs,cmap=plt.cm.RdBu_r,animated=True)
ugrid,newlons = shiftgrid(180.,u,longitudes,start=False)
vgrid,newlons = shiftgrid(180.,v,longitudes,start=False)
uproj,vproj,xx,yy = m.transform_vector(ugrid,vgrid,newlons,latitudes,31,31,returnxy=True,masked=True)
Q = m.quiver(xx,yy,uproj,vproj,scale=700)
qk = plt.quiverkey(Q, 0.1, 0.1, 20, '20 m/s', labelpos='W')
m.drawcoastlines(linewidth=1.5)
m.drawparallels(parallels)
m.drawmeridians(meridians)
cb = m.colorbar(CS2,"bottom", size="5%", pad="2%")
cb.set_label('hPa')
ax.set_title('SLP and Wind Vectors '+str(date.date()))
plt.show()
These have just been one time-step in a bigger dataset...
Use a local dataset for speed - Met Office CRUTEM4 (http://www.metoffice.gov.uk/hadobs/crutem4/data/gridded_fields/CRUTEM.4.4.0.0.anomalies.nc.gz)
First step is to turn a plot into a function call, then decorate it with a widget to scroll through the time-steps.
In [ ]:
from ipywidgets import interact, interactive, fixed
import ipywidgets as widgets
%matplotlib inline
In [ ]:
nc = Dataset('Data/CRUTEM.4.4.0.0.anomalies.nc')
lats = nc.variables['latitude'][:]
lons = nc.variables['longitude'][:]
lons, lats = np.meshgrid(lons,lats)
tind = nc.variables['time'][:]
@interact(index=(0, len(tind)))
def ftest(index):
basetime = datetime(1850, 1, 1)
date = basetime + timedelta(days=int(tind[index]))
tanom = nc.variables['temperature_anomaly'][index,:].squeeze()
fig = plt.figure(figsize=(10,10))
ax = fig.add_axes([0.05,0.05,0.9,0.9])
m = Basemap(projection='moll',lon_0=0,resolution='c')
m.drawcoastlines()
im1 = m.pcolormesh(lons, lats, tanom, shading='flat', cmap=cm.RdBu_r, latlon=True, vmin=-10, vmax=10)
m.drawparallels(np.arange(-90., 99., 30.))
m.drawmeridians(np.arange(-180., 180., 60.))
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
ax.set_title('{0} CRUTEM4 anom. (°C)'.format(date.date()))
plt.show()
In [ ]: