# 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 io import StringIO
from retrying import retry
%matplotlib inline
def princax(w):
# [theta,maj,min,wr] = princax(w)
#
# Input: w = complex vector time series (u+i*v)
#
# Output: theta = angle of maximum variance, (east == 0, north=90)
# maj = major axis of principal ellipse
# min = minor axis of principal ellipse
# wr = rotated time series, where real(wr) is aligned
# with the major axis.
cv = np.cov(np.real(w[:]),np.imag(w[:]))
#---------------------------------------------------------------------
# Find direction of maximum variance
#---------------------------------------------------------------------
theta = 0.5*np.arctan2(2*cv[1,0],(cv[0,0]-cv[1,1]))
#---------------------------------------------------------------------
# Find major and minor axis amplitudes
#---------------------------------------------------------------------
term1 = (cv[0,0]+cv[1,1])
term2 = np.sqrt((cv[0,0]-cv[1,1])**2 + 4*cv[1,0]**2)
maj = np.sqrt(.5*(term1+term2))
min = np.sqrt(.5*(term1-term2))
#---------------------------------------------------------------------
# Rotate into principal ellipse orientation
#---------------------------------------------------------------------
wr = w*np.exp(-1j*theta)
theta = theta*180./np.pi
return theta,maj,min,wr
# Specify which forecast model
#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
Massimo NW, 41.1405179, -69.0941178786
Massimo NE, 41.1405179863889, -68.7402883311111
Massimo SE, 40.8934875841667,-68.7402883311111
Massimo SW, 40.8934875841667,-69.0941178786111
'''
# Create a Pandas DataFrame from the station list
obs=pd.read_csv(StringIO(x.strip()), sep=",\s*",index_col='Station', engine='python')
obs
Lat | Lon | |
---|---|---|
Station | ||
Massimo NW | 41.140518 | -69.094118 |
Massimo NE | 41.140518 | -68.740288 |
Massimo SE | 40.893488 | -68.740288 |
Massimo SW | 40.893488 | -69.094118 |
# 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
@retry(stop_max_attempt_number=5, wait_fixed=3000)
def get_nc(url):
nc=netCDF4.Dataset(url)
return nc
nc = get_nc(url)
# 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 | |||
Massimo NW | 41.140518 | -69.094118 | 42427 |
Massimo NE | 41.140518 | -68.740288 | 43458 |
Massimo SE | 40.893488 | -68.740288 | 41829 |
Massimo SW | 40.893488 | -69.094118 | 38774 |
# get time values and convert to datetime objects
times = nc['time']
jd = netCDF4.num2date(times[:],times.units)
# get data for all time steps
nsta = len(obs)
u = np.ones((len(jd),nsta))
v = np.ones((len(jd),nsta))
s = np.ones((len(jd),nsta))
wr = np.ones((len(jd),nsta))
stheta = np.ones(nsta)
smaj = np.ones(nsta)
smin = np.ones(nsta)
#klev = 0 # surface layer
klev = -1 # bottom layer in FVCOM is 1 m above bottom
for i in range(nsta):
u[:,i] = nc['u'][:,klev,obs['0-Based Index'][i]]
v[:,i] = nc['v'][:,klev,obs['0-Based Index'][i]]
s[:,i] = np.sqrt(u[:,i]**2 + v[:,i]**2)
[stheta[i], smaj[i], smin[i], wr[:,i]] = princax(u[:,i] + 1j * v[:,i])
C:\Users\rsignell\miniconda3\envs\IOOS\lib\site-packages\ipykernel\__main__.py:18: ComplexWarning: Casting complex values to real discards the imaginary part
# make a DataFrame out of the interpolated time series at each location
zvals=pd.DataFrame(v,index=jd,columns=obs.index)
# list out a few values
zvals.head()
Station | Massimo NW | Massimo NE | Massimo SE | Massimo SW |
---|---|---|---|---|
2018-06-03 00:00:00.000 | -0.084284 | -0.113184 | -0.147742 | -0.179091 |
2018-06-03 01:01:52.500 | 0.097977 | 0.080652 | 0.047547 | 0.006214 |
2018-06-03 01:58:07.500 | 0.229463 | 0.223774 | 0.193057 | 0.161910 |
2018-06-03 03:00:00.000 | 0.302423 | 0.305330 | 0.265874 | 0.234719 |
2018-06-03 04:01:52.500 | 0.278393 | 0.305741 | 0.261742 | 0.233263 |
# 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['u'].units
# Add new column to DataFrame
b[zmax_heading]=zvals.max()
b
time of max water level (UTC) | zmax (meters s-1) | |
---|---|---|
Station | ||
Massimo NW | 2018-06-08 19:58:07.500 | 0.341587 |
Massimo NE | 2018-06-08 19:58:07.500 | 0.346320 |
Massimo SE | 2018-06-08 19:58:07.500 | 0.313111 |
Massimo SW | 2018-06-08 19:58:07.500 | 0.265076 |
# plotting at DataFrame is easy!
ax=zvals.plot(figsize=(16,4),grid=True,title=('NECOFS bottom northward current 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 DataFrame out of the interpolated time series at each location
zvals=pd.DataFrame(u,index=jd,columns=obs.index)
# plotting at DataFrame is easy!
ax=zvals.plot(figsize=(16,4),grid=True,title=('NECOFS eastward 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 DataFrame out of the interpolated time series at each location
svals=pd.DataFrame(s,index=jd,columns=obs.index)
# plotting at DataFrame is easy!
ax=svals.plot(figsize=(16,4),grid=True,title=('NECOFS bottom current 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));
print(u.mean(),v.mean())
0.005482110726006795 -0.07266947751431244
stheta
array([ 78.03716566, -87.02170566, -89.35171984, 78.43577375])