# Plot forecast water levels from NECOFS model from list of lon,lat locations
# (uses the nearest point, no interpolation)
import netCDF4
import datetime as dt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from StringIO import StringIO
%matplotlib inline
#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'
def dms2dd(d,m,s):
return d+(m+s/60.)/60.
In http://openscholar.purchase.edu/george_kraemer/files/2015_kim_et_al._kelp.pdf there are three kelp sites are listed as:
But since BRE site is obviously wrong (80 and 87 minutes > 60 minutes), we used https://openscholar.purchase.edu/george_kraemer/files/2014_kim_et_al._grac_aquaculture.pdf where BRE is listed as
print(dms2dd(40,48.047,0))
print(-dms2dd(73,52.164,0))
40.8007833333 -73.8694
x = '''
Station, Lat, Lon
Fairfield, 41.11470, -73.25461
Branford, 41.21866, -73.95117
Bronx River Estuary, 40.80078, -73.8694
'''
# Create a Pandas DataFrame
obs=pd.read_csv(StringIO(x.strip()), sep=",\s*",index_col='Station')
-c:2: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support regex separators (separators > 1 char and different from '\s+' are interpreted as regex); you can avoid this warning by specifying engine='python'.
obs
Lat | Lon | |
---|---|---|
Station | ||
Fairfield | 41.11470 | -73.25461 |
Branford | 41.21866 | -73.95117 |
Bronx River Estuary | 40.80078 | -73.86940 |
# find the indices of the points in (x,y) closest to the points in (xi,yi)
def nearxy(x,y,xi,yi):
ind = np.ones(len(xi),dtype=int)
for i in np.arange(len(xi)):
dist = np.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['lonc'][:],nc['latc'][:],obs['Lon'],obs['Lat'])
obs
Lat | Lon | 0-Based Index | |
---|---|---|---|
Station | |||
Fairfield | 41.11470 | -73.25461 | 76371 |
Branford | 41.21866 | -73.95117 | 97584 |
Bronx River Estuary | 40.80078 | -73.86940 | 47635 |
# 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
klev = 0 #FVCOM layer: 0 = surface, -1 = bottom
nsta = len(obs)
z = np.ones((len(jd),nsta))
for i in range(nsta):
u = nc['u'][:,klev,obs['0-Based Index'][i]]
v = nc['v'][:,klev,obs['0-Based Index'][i]]
z[:,i] = np.sqrt(u**2 + v**2)
# 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 | Fairfield | Branford | Bronx River Estuary |
---|---|---|---|
2017-04-11 00:00:00.000 | 0.324626 | 0.050522 | 1.012293 |
2017-04-11 01:01:52.500 | 0.361795 | 0.074158 | 1.049958 |
2017-04-11 01:58:07.500 | 0.273000 | 0.129365 | 0.799587 |
2017-04-11 03:00:00.000 | 0.107714 | 0.145893 | 0.175020 |
2017-04-11 04:01:52.500 | 0.081222 | 0.088899 | 0.773422 |
# plotting at DataFrame is easy!
ax=zvals.plot(figsize=(16,4),grid=True,title=('NECOFS Forecast Water Speed from %s Forecast' % model),legend=False);
# read units from dataset for ylabel
plt.ylabel(nc['u'].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));
# make a new DataFrame of maximum water speed at all stations
b=pd.DataFrame(zvals.idxmax(),columns=['time of max water speed (UTC)'])
# create heading for new column containing max water speed
zmax_heading='umax (%s)' % nc['u'].units
# Add new column to DataFrame
b[zmax_heading]=zvals.max()
b
time of max water speed (UTC) | umax (meters s-1) | |
---|---|---|
Station | ||
Fairfield | 2017-04-12 13:58:07.500 | 0.404602 |
Branford | 2017-04-13 10:01:52.500 | 0.270949 |
Bronx River Estuary | 2017-04-12 06:00:00.000 | 1.169032 |