# 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.
dms2dd(41,30.2,0)
41.50333333333333
-dms2dd(71,19.6,0)
-71.32666666666667
x = '''
Station, Lat, Lon
Falmouth Harbor, 41.541575, -70.608020
Sage Lot Pond, 41.554361, -70.505611
'''
x = '''
Station, Lat, Lon
Boston, 42.368186, -71.047984
Carolyn Seep Spot, 39.8083, -69.5917
Falmouth Harbor, 41.541575, -70.608020
'''
# 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
Newport, 41.50333, -71.3266
'''
# 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 | ||
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 |
Newport | 41.503330 | -71.326600 |
# 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['lon'][:],nc['lat'][:],obs['Lon'],obs['Lat'])
obs
Lat | Lon | 0-Based Index | |
---|---|---|---|
Station | |||
Boston | 42.368186 | -71.047984 | 51703 |
Scituate Harbor | 42.199447 | -70.720090 | 47715 |
Scituate Beach | 42.209973 | -70.724523 | 47710 |
Falmouth Harbor | 41.541575 | -70.608020 | 44738 |
Marion | 41.689008 | -70.746576 | 49963 |
Marshfield | 42.108480 | -70.648691 | 43039 |
Provincetown | 42.042745 | -70.171180 | 41278 |
Sandwich | 41.767990 | -70.466219 | 46662 |
Hampton Bay | 42.900103 | -70.818510 | 48804 |
Gloucester | 42.610253 | -70.660570 | 39565 |
Newport | 41.503330 | -71.326600 | 43767 |
# 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 = np.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 | Newport |
---|---|---|---|---|---|---|---|---|---|---|---|
2017-03-11 00:00:00.000 | 0.253245 | 0.372455 | 0.384904 | 0.297232 | 0.871409 | 0.377139 | 0.317074 | 0.370238 | 0.465665 | 0.423959 | 0.855790 |
2017-03-11 01:01:52.500 | 1.045947 | 1.024449 | 1.046528 | 0.137273 | 0.449498 | 1.027487 | 0.986655 | 1.028728 | 1.128143 | 1.064781 | 0.554442 |
2017-03-11 01:58:07.500 | 1.474221 | 1.494632 | 1.516271 | -0.066920 | -0.070274 | 1.516882 | 1.535385 | 1.576048 | 1.530783 | 1.490083 | 0.081553 |
2017-03-11 03:00:00.000 | 1.687645 | 1.624390 | 1.645241 | -0.137259 | -0.499917 | 1.652100 | 1.734829 | 1.758991 | 1.595022 | 1.583244 | -0.353817 |
2017-03-11 04:01:52.500 | 1.368101 | 1.307220 | 1.323382 | -0.093599 | -0.709073 | 1.341266 | 1.466981 | 1.456549 | 1.241830 | 1.246800 | -0.670149 |
# 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
plt.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['Boston'].max()
1.9948580265045166
# plotting at DataFrame is easy!
ax=zvals['Newport'].plot(figsize=(16,4),grid=True,title=('NECOFS Forecast Water Level from %s Forecast' % model),legend=False);
# read units from dataset for ylabel
plt.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));
# make a new DataFrame of maximum water levels at all stations
b=pd.DataFrame(zvals.idxmax(),columns=['time of max water level (UTC)'])
# create heading for new column containing max water level
zmax_heading='zmax (%s)' % nc['zeta'].units
# Add new column to DataFrame
b[zmax_heading]=zvals.max()
b
time of max water level (UTC) | zmax (meters) | |
---|---|---|
Station | ||
Boston | 2017-03-14 18:00:00.000 | 1.994858 |
Scituate Harbor | 2017-03-14 18:00:00.000 | 1.787220 |
Scituate Beach | 2017-03-11 15:00:00.000 | 1.765060 |
Falmouth Harbor | 2017-03-14 15:00:00.000 | 0.648843 |
Marion | 2017-03-15 03:00:00.000 | 1.359447 |
Marshfield | 2017-03-11 15:00:00.000 | 1.776518 |
Provincetown | 2017-03-11 15:00:00.000 | 1.837048 |
Sandwich | 2017-03-11 15:00:00.000 | 1.868223 |
Hampton Bay | 2017-03-14 16:58:07.500 | 1.765515 |
Gloucester | 2017-03-14 16:58:07.500 | 1.716410 |
Newport | 2017-03-14 15:00:00.000 | 1.077989 |