Contour SST from OPeNDAP server using Python

I often want to make a quick plot of the observed sea surface temperature for a specified region to compare with output from ocean circulation models. I also often want the data to be of relatively high-resolution so that you can see the details and structures along coastlines as well as the details of gyres and meanders in the ocean. This usually means that you either need to download a large dataset to your local storage or you can use OPeNDAP (Open-source Project for a Network Data Access Protocol). When you combine OPeNDAP with Python you can easily extract the data you need for a given region and time-period and do the plotting quite easily using Matplotlib. Below is a gallery of the images I created using the approach described on this page.

SST plot created by getCortad.py

The following example will extract SST data from an online archive called CoRTAD available from NODC NOAA. The archive includes blended night/day Pathfinder 5 SST data that have been interpolated to fill gaps caused by clouds etc. The dataset provides SST at 4×4 km resolution globally and is a fantastic source for ocean temperature starting in 1982 and ending in 2009. If you use this dataset make sure to reference NOAA NODC and also a peer-review article that describes the dataset: Selig et al. 2010 (787)New insights into global patterns of ocean temperature anomalies: implications for coral reef health and management. Global Ecology and Biogeography 19, 397–411.

For this example we will extract SST for the North Sea region. The CoRTAD data are provided as three dimensional cubes of data that give information of SST (and other variables) for a specified geographial region (a tile), and the third dimension represents time. If you look at the tile-map of how the CoRTAD data are organized (click here) you will notice that the North Sea region is represented by two separate tiles (1.7 and 1.8). This means that we need to extract both tiles and then concentrate them afterwards to create a single plot of SST for a given date. This is not straightforward, but if you remember a couple of steps its alright. First, get the dimensions of the different tiles by extracting the longitude and latitudes:

base="http://data.nodc.noaa.gov/opendap/cortad/Version3/"
file1="cortadv3_row01_col07.h5"
filename1=base+file1
cdf1=Dataset(filename1)
longitude1=np.squeeze(cdf1.variables["Longitude"][:])
latitude1=np.squeeze(cdf1.variables["Latitude"][:])

Repeat this for tile number two and then get the time variable

time=np.squeeze(cdf1.variables["Time"][:])

Time in the CoRTAD file is represented as “days since 1982-01-01-00:00:00”. Python provides a very useful module

 import datetime

which can help you modify the array of days since 1982 (time array). Define a reference date

 refDate = datetime.datetime(1982,1,1,0,0,0)

and then access the dates in the time array by using

for t in range(len(time)):
currentTime=refDate + datetime.timedelta(days=time[t])

This gives you the ability to loop through the time array and find the index representing the date you want to extract. Once you have the index, then extract the data for that time-index for both of the tiles that represent the North Sea region. As this is HDF data we need to account for offset in the variables and finally we need to concentate the data for both regions into one array that is useful for plotting.

filledSST1=cdf1.variables["FilledSST"][:,:,myIndex]
filledSST2=cdf2.variables["FilledSST"][:,:,myIndex]
offset1=cdf1.variables["FilledSST"].__getattribute__('add_offset')
offset2=cdf2.variables["FilledSST"].__getattribute__('add_offset')
filledSST1=filledSST1 + abs(offset1)
filledSST1=np.rot90(filledSST1,3) # rotate 90 degrees 3 times
filledSST1=np.fliplr(filledSST1) # then flip left right
filledSST2=filledSST2 + abs(offset2)
filledSST2=np.rot90(filledSST2,3) # rotate 90 degrees 3 times
filledSST2=np.fliplr(filledSST2) # then flip left right
filledSST=concatenate((filledSST1,filledSST2),axis=direction)

We also need to add the longitude and latitude arrays together

if direction==1:
longitude=concatenate((longitude1,longitude2),axis=direction)
latitude=latitude1
else:
longitude=longitude1
latitude=concatenate((latitude1,latitude2),axis=direction)

Once you have the data concentated you are ready to contour the results. I usually use the “contourf” function from Matplotlib as I think it looks great and is visually easy to understand. The main part of the drawing routine is:

levels = np.arange(-2.0, 20.0, 0.5)
CS1 = map.contourf(x,y,filledSST,levels,
cmap=cm.get_cmap('jet',len(levels)-1),
extend='both',
alpha=1.0,
origin='lower',
rasterized=True)

Here, levels represents the number of contours to draw for the area. I selected to draw every 0.5 degrees Celsius from -2.0C to 20C, which is quite representable for the North Sea during the year. The final result is plotted and saved to file. The full code for this example is downloadable below. One problem I had with this script was that it would require a lot of memory if you choose to plot a lot of years and months. If anyone has suggestions for how to avoid the extreme memory usage or have ideas for how to optimize this script I would appreciate comments or an email.

download getCortad.py (842)

import os, sys, datetime, string
import numpy as np
from netCDF4 import Dataset
import numpy.ma as ma
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, shiftgrid, NetCDFFile
from pylab import *

__author__ = 'Trond Kristiansen'
__email__ = 'trond.kristiansen(at)imr.no'
__created__ = datetime.datetime(2011, 6, 25)
__modified__ = datetime.datetime(2012, 4, 10)
__version__ = "1.0"
__status__ = "Development, 25.6.2011, 10.04.2012"

def makeMap(latStart,latEnd,lonStart,lonEnd,longitude,latitude,filledSST,name,dateString):

plt.figure(figsize=(12,10),frameon=False)

if lonStart< 0 and lonEnd < 0: lon_0= - (abs(lonEnd)+abs(lonStart))/2.0 elif lonStart > 0 and lonEnd > 0:
lon_0=(abs(lonEnd)+abs(lonStart))/2.0
else:
lon_0=((lonEnd)+(lonStart))/2.0

map = Basemap(llcrnrlat=latStart,urcrnrlat=latEnd,\
llcrnrlon=lonStart,urcrnrlon=lonEnd,\
rsphere=(6378137.00,6356752.3142),\
resolution='l',area_thresh=1000.,projection='lcc',\
lat_1=latStart,lon_0=lon_0)

lon,lat=meshgrid(longitude,latitude)

x, y = map(lon,lat)
map.drawcoastlines()
map.fillcontinents(color='grey')
map.drawcountries()
#map.bluemarble()

levels = np.arange(-2.0, 20.0, 0.5)

CS1 = map.contourf(x,y,filledSST,levels,
cmap=cm.get_cmap('jet',len(levels)-1),
extend='both',
alpha=1.0,
origin='lower',
rasterized=True)

CS1.axis='tight'
map.drawmeridians(np.arange(lon.min(),lon.max(),10),labels=[0,0,0,1])
map.drawparallels(np.arange(lat.min(),lat.max(),4),labels=[1,0,0,0])

plt.colorbar(CS1,orientation='vertical',extend='both', shrink=0.5)

plt.title('CoRTAD FilledSST - %s'%(dateString))
plotfile='figures/SST_'+str(dateString)+'_'+str(name)+'.png'
plt.savefig(plotfile)
plt.show()

def extratCORTAD(name,latStart,latEnd,lonStart,lonEnd):

"""
Info on the different tiles used to identoify a region is found here:
http://www.nodc.noaa.gov/SatelliteData/Cortad/TileMap.jpg
"""

base="http://data.nodc.noaa.gov/opendap/cortad/Version3/"

refDate = datetime.datetime(1982,1,1,0,0,0)

if name=="North Sea":
file1="cortadv3_row01_col07.h5"
file2="cortadv3_row01_col08.h5"
tiles=2
direction=1 # add the tiles

filename1=base+file1
filename2=base+file2
print "Extrating data from openDAP: %s"%(filename1)
print "Extrating data from openDAP: %s"%(filename2)

cdf1=Dataset(filename1)
cdf2=Dataset(filename2)

longitude1=np.squeeze(cdf1.variables["Longitude"][:])
latitude1=np.squeeze(cdf1.variables["Latitude"][:])

longitude2=np.squeeze(cdf2.variables["Longitude"][:])
latitude2=np.squeeze(cdf2.variables["Latitude"][:])
time=np.squeeze(cdf2.variables["Time"][:])
print "Finished extracting the grid"

wantedYears=[1990]
wantedMonths=[1,2,3,4,5,6,7,8,9,10,11,12]
wantedMonths=[1]

for t in range(len(time)):
currentTime=refDate + datetime.timedelta(days=time[t])
#print currentTime
if currentTime.year in wantedYears and currentTime.month in wantedMonths:
print "Current time inspecting: %s"%(currentTime)
myIndex=t; dateString="%s"%(currentTime)

filledSST1=cdf1.variables["FilledSST"][:,:,myIndex]
filledSST2=cdf2.variables["FilledSST"][:,:,myIndex]

offset1=cdf1.variables["FilledSST"].__getattribute__('add_offset')
offset2=cdf2.variables["FilledSST"].__getattribute__('add_offset')

filledSST1=filledSST1 + abs(offset1)
filledSST1=np.rot90(filledSST1,3) # rotate 90 degrees 3 times
filledSST1=np.fliplr(filledSST1) # then flip left right

filledSST2=filledSST2 + abs(offset2)
filledSST2=np.rot90(filledSST2,3) # rotate 90 degrees 3 times
filledSST2=np.fliplr(filledSST2) # then flip left right

filledSST=concatenate((filledSST1,filledSST2),axis=direction)

if direction==1:
longitude=concatenate((longitude1,longitude2),axis=direction)
latitude=latitude1
else:
longitude=longitude1
latitude=concatenate((latitude1,latitude2),axis=direction)

print "Created matrix for SST for area: %s"%(name)
print "Matrix size: ",filledSST.shape
print "Longitude size: ",longitude.shape
print "Latitude size: ",latitude.shape
print "Long min: %s Long max: %s"%(longitude.min(),longitude.max())
print "Lat min: %s Lat max: %s"%(latitude.min(),latitude.max())
print "------------------------------\n"
makeMap(latStart,latEnd,lonStart,lonEnd,longitude,latitude,filledSST,name,dateString)

def main():
latStart = 43
latEnd = 68
lonStart =-19
lonEnd = 33
name='North Sea'
extratCORTAD(name,latStart,latEnd,lonStart,lonEnd)

if __name__ == "__main__":
main()

2 thoughts on “Contour SST from OPeNDAP server using Python”

  1. I appreciated for your valuable script. Thank you very much.
    But I would like to ask you about my problem when I run it shown below: from mpl_toolkits.basemap import Basemap, shiftgrid, NetCDFFile ImportError: cannot import name NetCDFFile

    Could you help me ?

    best regard

    dadang s

  2. Hi. Just remove importing NetCDFFile (mpl_toolkits.basemap import Basemap, shiftgrid) and it should work. I used to use NetCDFFile but now I use the NetCDF4 module instead. Cheers, Trond

Leave a Reply

Your email address will not be published. Required fields are marked *