Accessing oceans of model data: Oceanographic analysis using IPython Notebook & ArcGIS Python interface

Rich Signell, Research Oceanographer

US Geological Survey, Woods Hole, MA

[email protected]

http://bitly.com/bundles/rsignell/1


The Problem: The data is big. This 33 year ocean hindcast dataset is 33TB. It is also unstructured (on a triangular mesh). How can researchers and environmental managers access this data?


Step 1. Make the NetCDF output files accessible via OPeNDAP (Most commonly used Web Service for climate, atmospheric and ocean model output

In [1]:
from IPython.core.display import HTML
HTML('<iframe src=http://www.smast.umassd.edu:8080/thredds/hindcasts.html?dataset=fvcom/hindcasts/30yr_gom3 width=1000 height=600></iframe>')
Out[1]:

Step 2. Use NetCDF4-Python in Ipython Notebook for quick visualization and analysis, working close to the data

In [2]:
import netCDF4
import matplotlib.pyplot as plt
import matplotlib.tri as tri
In [3]:
url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3/mean'
In [4]:
nc = netCDF4.Dataset(url)
x = nc.variables['lon'][:]
y = nc.variables['lat'][:]
var = nc.variables['temp'][-1,-1,:]   # [time, level, node]
# read connectivity array
nv= nc.variables['nv'][:,:]
nv= nv.T - 1  # 1-based => 0-based indexing of nodes
nc.close()
In [5]:
# define TRI object
triang = tri.Triangulation(x, y, triangles=nv)
In [9]:
# Tricontourf is fast 
plt.figure(figsize=(18,10))
plt.gca().set_aspect(1./cos(42*pi/180))
plt.tricontourf(triang, var,levels=arange(0,20,1));
plt.colorbar();
#plt.tricontour(triang, var, colors='k',levels=arange(0,20,1))
plt.title('Bottom Temperature, December, 2010');
In [7]:
# Tripcolor is slow
plt.figure(figsize=(18,10))
plt.gca().set_aspect('equal')
plt.tripcolor(x, y, nv, var, shading='faceted')
plt.axis([-71.2, -70.5, 42.0, 42.5])
plt.colorbar()
plt.title('Bottom Temperature in December, 2010')
plt.xlabel('Longitude (degrees)')
plt.ylabel('Latitude (degrees)')
Out[7]:
<matplotlib.text.Text at 0x3c36250>

Step 3. Use ArcPy in ArcGIS 10.1 with NetCDF4-Python to access and perform GIS analysis on the desktop

In [8]:
from IPython.core.display import Image
Image('http://epi.whoi.edu/esr/Picture_ESRI.png')
Out[8]: