for Chris Ackland - sst ramssa plot

ramssa.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Author: Laurent Besnard
# Institute: IMOS / eMII
# email address: laurent.besnard@utas.edu.au

from netCDF4 import Dataset, num2date
from datetime import date
from numpy import meshgrid
from matplotlib.pyplot import (figure, pcolor, get_cmap, colorbar,
                                       xlabel, ylabel, title, show)

srsUrl         = 'http://opendap.jpl.nasa.gov:80/opendap/allData/ghrsst/data/L4/AUS/ABOM/RAMSSA_09km/2016/001/20160101-ABOM-L4HRfnd-AUS-v01-fv01_0-RAMSSA_09km.nc.bz2'
geoBoundaryBox = [165, 181, -50 ,-30] # [lon min , lonmax , latmin , latmax] . New Zealand
srsDataObj     = Dataset(srsUrl)

# load all the latitude and longitude values to find the indexes which will match geoBoundaryBox in order to subset the variables
sst  = srsDataObj.variables['analysed_sst']
lon  = srsDataObj.variables['lon']
lat  = srsDataObj.variables['lat']
mask = srsDataObj.variables['mask']
time = srsDataObj.variables['time']

latValues              = lat[:]
lonValues              = lon[:]
lonValues[lonValues<0] = lonValues[lonValues<0] + 360 # modify the longitude values which are across the 180th meridian
timeValues             = num2date(time[:], time.units)

indexLatSubSelect   = (latValues < geoBoundaryBox[3]) & (latValues > geoBoundaryBox[2])
indexLonSubSelect   = (lonValues < geoBoundaryBox[1]) & (lonValues > geoBoundaryBox[0])

sstValuesSubSelect  = sst[0, indexLatSubSelect,indexLonSubSelect]
latValuesSubSelect  = latValues[indexLatSubSelect]
lonValuesSubSelect  = lonValues[indexLonSubSelect]
maskValuesSubSelect = mask[0, indexLatSubSelect, indexLonSubSelect]  # land mask

[lon_mesh, lat_mesh] = meshgrid(lonValuesSubSelect, latValuesSubSelect) # we create a matrix of similar size to be used afterwards with pcolor

figure1 =  figure(num=None, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')
pcolor(lon_mesh, lat_mesh , maskValuesSubSelect == 2, cmap=get_cmap('gray'))

pcolor(lon_mesh, lat_mesh ,sstValuesSubSelect)
title('%s - %s' % (srsDataObj.title, timeValues[0].strftime('%Y-%m-%d %H:%M:%S')))
xlabel('%s in %s' % (lon.long_name, lon.units))
ylabel('%s in %s' % (lat.long_name, lat.units))

cbar = colorbar()
cbar.ax.set_ylabel('%s \n in %s' % (sst.long_name, sst.units))
show()