# 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
#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'
# 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
'''
# Create a Pandas DataFrame
obs=pd.read_csv(StringIO(x.strip()), sep=",\s*",index_col='Station')
obs
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 |
# 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
# open NECOFS remote OPeNDAP dataset
nc=netCDF4.Dataset(url).variables
# find closest NECOFS nodes to station locations
obs['0-Based Index'] = nearxy(nc['lon'][:],nc['lat'][:],obs['Lon'],obs['Lat'])
obs
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 |
# get time values and convert to datetime objects
times = nc['time']
jd = netCDF4.num2date(times[:],times.units)
# 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]]
# make a DataFrame out of the interpolated time series at each location
zvals=pd.DataFrame(z,index=jd,columns=obs.index)
# list out a few values
zvals.head()
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 |
# 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));
# what is the maximum water level at Scituate over this period?
zvals['Scituate Harbor'].max()
1.8255265951156616
# 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()
b
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 |