Extract NECOFS water levels using NetCDF4-Python and analyze/visualize with Pandas

In [1]:
# Plot forecast water levels from NECOFS model from list of lon,lat locations
# (uses the nearest point, no interpolation)
from pylab import *
import netCDF4
import datetime as dt
import pandas as pd
from StringIO import StringIO
In [2]:
#NECOFS MassBay grid
model='Massbay'
url='http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc'
# GOM3 Grid
#model='GOM3'
#url='http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc'
In [3]:
# Enter desired (Station, Lat, Lon) values here:
x = '''
Station, Lat, Lon
Boston,             42.368186, -71.047984
Scituate Harbor,    42.199447, -70.720090
Scituate Beach,     42.209973, -70.724523
Falmouth Harbor,    41.541575, -70.608020
Marion,             41.689008, -70.746576
Marshfield,         42.108480, -70.648691
Provincetown,       42.042745, -70.171180
Sandwich,           41.767990, -70.466219
Hampton Bay,        42.900103, -70.818510
Gloucester,         42.610253, -70.660570
'''
In [4]:
# Create a Pandas DataFrame
obs=pd.read_csv(StringIO(x.strip()), sep=",\s*",index_col='Station')
In [5]:
obs
Out[5]:
Lat Lon
Station
Boston 42.368186 -71.047984
Scituate Harbor 42.199447 -70.720090
Scituate Beach 42.209973 -70.724523
Falmouth Harbor 41.541575 -70.608020
Marion 41.689008 -70.746576
Marshfield 42.108480 -70.648691
Provincetown 42.042745 -70.171180
Sandwich 41.767990 -70.466219
Hampton Bay 42.900103 -70.818510
Gloucester 42.610253 -70.660570
In [6]:
# find the indices of the points in (x,y) closest to the points in (xi,yi)
def nearxy(x,y,xi,yi):
    ind=ones(len(xi),dtype=int)
    for i in arange(len(xi)):
        dist=sqrt((x-xi[i])**2+(y-yi[i])**2)
        ind[i]=dist.argmin()
    return ind
In [7]:
# open NECOFS remote OPeNDAP dataset 
nc=netCDF4.Dataset(url).variables
In [8]:
# find closest NECOFS nodes to station locations
obs['0-Based Index'] = nearxy(nc['lon'][:],nc['lat'][:],obs['Lon'],obs['Lat'])
obs
Out[8]:
Lat Lon 0-Based Index
Station
Boston 42.368186 -71.047984 90913
Scituate Harbor 42.199447 -70.720090 37964
Scituate Beach 42.209973 -70.724523 28474
Falmouth Harbor 41.541575 -70.608020 47470
Marion 41.689008 -70.746576 49654
Marshfield 42.108480 -70.648691 24272
Provincetown 42.042745 -70.171180 26595
Sandwich 41.767990 -70.466219 38036
Hampton Bay 42.900103 -70.818510 13022
Gloucester 42.610253 -70.660570 22082
In [9]:
# get time values and convert to datetime objects
times = nc['time']
jd = netCDF4.num2date(times[:],times.units)
In [10]:
# get all time steps of water level from each station
nsta=len(obs)
z=ones((len(jd),nsta))
for i in range(nsta):
    z[:,i] = nc['zeta'][:,obs['0-Based Index'][i]]
    
In [11]:
# make a DataFrame out of the interpolated time series at each location
zvals=pd.DataFrame(z,index=jd,columns=obs.index)
In [12]:
# list out a few values
zvals.head()
Out[12]:
Station Boston Scituate Harbor Scituate Beach Falmouth Harbor Marion Marshfield Provincetown Sandwich Hampton Bay Gloucester
2013-04-11 00:00:00 -0.445120 0.486000 -0.429675 0.509572 1.056936 -0.443456 -0.528928 -0.503050 -0.289000 -0.413870
2013-04-11 01:01:53 0.129476 0.486000 0.091739 0.587705 1.235945 0.145963 0.168091 0.234113 0.128534 0.182251
2013-04-11 01:58:08 0.799068 0.823351 0.849990 0.551944 0.956099 0.820433 0.758173 0.811803 0.922069 0.878746
2013-04-11 03:00:00 1.572014 1.410582 1.421661 0.419279 0.584492 1.410825 1.370721 1.419978 1.426254 1.435683
2013-04-11 04:01:53 1.790708 1.766814 1.801998 0.291546 0.124641 1.798690 1.825977 1.861521 1.769746 1.796204
In [13]:
# plotting at DataFrame is easy!
ax=zvals.plot(figsize=(16,4),grid=True,title=('NECOFS Forecast Water Level from %s Forecast' % model),legend=False);
# read units from dataset for ylabel
ylabel(nc['zeta'].units)
# plotting the legend outside the axis is a bit tricky
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5));
In [14]:
# what is the maximum water level at Scituate over this period?
zvals['Scituate Harbor'].max()
Out[14]:
1.8255265951156616
In [15]:
# make a new DataFrame of maximum water levels at all stations
b=pd.DataFrame(zvals.idxmax(),columns=['time of max water level (UTC)'])
b['zmax (%s)' % nc['zeta'].units]=zvals.max()
In [16]:
b
Out[16]:
time of max water level (UTC) zmax (meters)
Station
Boston 2013-04-11 04:58:08 1.938324
Scituate Harbor 2013-04-12 04:58:08 1.825527
Scituate Beach 2013-04-12 04:58:08 1.823399
Falmouth Harbor 2013-04-12 01:58:08 0.659575
Marion 2013-04-11 01:01:53 1.235945
Marshfield 2013-04-12 04:58:08 1.831563
Provincetown 2013-04-12 04:58:08 1.866504
Sandwich 2013-04-12 04:58:08 1.901897
Hampton Bay 2013-04-12 04:58:08 1.816520
Gloucester 2013-04-12 04:58:08 1.797317