# 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
import time
# currents
urls={}
urls['FVCOM Currents']='http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3'
urls['FVCOM Waves']='http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/wave_gom3'
sta = {}
sta['Outer Fall'] = (43.30, -68.65)
sta['Great South Channel'] = (41.29, -69.18)
nsta=len(sta)
loni=[]
lati=[]
names=[]
for ll in sta.itervalues():
lati.append(ll[0])
loni.append(ll[1])
for station_name, ll in sta.iteritems():
print '%s %9.4f, %9.4f' % (station_name,ll[0],ll[1])
names.append(station_name)
Great South Channel 41.2900, -69.1800 Outer Fall 43.3000, -68.6500
# 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(loni),dtype=int)
for i in arange(len(xi)):
dist=sqrt((x-xi[i])**2+(y-yi[i])**2)
ind[i]=dist.argmin()
return ind
nc=netCDF4.Dataset(urls['FVCOM Waves']);
# don't find closest cells if they already exist.
# so if you change the stations, you need to restart the kernel or do "del ind"
if 'lon' not in locals():
lon=nc.variables['lon'][:]
lat=nc.variables['lat'][:]
ind = nearxy(lon,lat,loni,lati)
ind
array([21937, 25660])
# get the time values and convert to datetime objects
times = nc.variables['time']
jdwave = netCDF4.num2date(times[:],times.units)
# get all time steps of waves from each station
hs=ones((len(jdwave),len(loni)))
for i in range(len(loni)):
hs[:,i] = nc.variables['hs'][:,ind[i]]
if 'lonc' not in locals():
lonc=nc.variables['lonc'][:]
latc=nc.variables['latc'][:]
ind = nearxy(lonc,latc,loni,lati)
ind
array([43629, 50361])
# get all time steps of waves from each station
uwind=ones((len(jdwave),len(loni)))
vwind=ones((len(jdwave),len(loni)))
time1=time.time()
for i in range(len(loni)):
uwind[:,i] = nc.variables['uwind_speed'][:,ind[i]]
vwind[:,i] = nc.variables['vwind_speed'][:,ind[i]]
time2=time.time()
print 'elapsed seconds = %f' % (time2-time1)
elapsed seconds = 284.227118
# make the time series plot, with nicely formatted labels
MyDateFormatter = DateFormatter('%b %d\n %Y')
fig = plt.figure(figsize=(22,6), dpi=80)
ax1 = fig.add_subplot(111)
#jd_est = jd - dt.timedelta(hours=5)
# plot entire hindcast/forecast
ax1.plot(jdwave,hs) # convert from meters to feet
# plot only more recent levels (from step 72 on...)
#ax1.plot(jd_est[72:],z[72:,:]/.3048) # convert from meters to feet
ax1.xaxis.set_major_formatter(MyDateFormatter)
ax1.grid(True)
setp(gca().get_xticklabels(), rotation=0, horizontalalignment='center')
#plt.title('NECOFS Forecast Waves from %s Forecast' % model)
ax1.set_ylabel('m')
ax1.legend(names,loc='upper left')
<matplotlib.legend.Legend at 0x3990a10>
wind_spd=sqrt(uwind*uwind+vwind*vwind)
# make the time series plot, with nicely formatted labels
MyDateFormatter = DateFormatter('%b %d\n %Y')
fig = plt.figure(figsize=(22,6), dpi=80)
ax1 = fig.add_subplot(111)
#jd_est = jd - dt.timedelta(hours=5)
# plot entire hindcast/forecast
ax1.plot(jdwave,wind_spd) # convert from meters to feet
# plot only more recent levels (from step 72 on...)
#ax1.plot(jd_est[72:],z[72:,:]/.3048) # convert from meters to feet
ax1.xaxis.set_major_formatter(MyDateFormatter)
ax1.grid(True)
setp(gca().get_xticklabels(), rotation=0, horizontalalignment='center')
#plt.title('NECOFS Forecast Waves from %s Forecast' % model)
ax1.set_ylabel('m/s')
ax1.legend(names,loc='upper left')
<matplotlib.legend.Legend at 0x3081fd0>
data = {}
for i in range(len(sta)):
data[sta.keys()[i]]=wind_spd[:,i]
df = pd.DataFrame(data,index=jdwave)
fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(111)
df.plot(ax=ax)
#df.plot(ax=ax,title='%s (Lon=%.4f, Lat=%.4f)' % (h.standard_name, lon, lat))
ax.set_ylabel(nc.variables['uwind_speed'].units)
<matplotlib.text.Text at 0x305ce10>
nco = netCDF4.Dataset('/usgs/data1/notebook/kemp_wind_and_wave.nc', 'w',clobber=True)
nco.createDimension('time', None)
nco.createDimension('station',nsta)
timeo = nco.createVariable('time', 'f4', ('time'))
lono = nco.createVariable('lon', 'f4', ('station'))
lato = nco.createVariable('lat', 'f4', ('station'))
hso = nco.createVariable('hsig', 'f4', ('time', 'station'))
uo = nco.createVariable('uwind', 'f4', ('time', 'station'))
vo = nco.createVariable('vwind', 'f4', ('time', 'station'))
nco.variables['time'].units = times.units
nco.variables['lon'].units='degrees_east'
nco.variables['lat'].units='degrees_north'
nco.variables['hsig'].units='m'
nco.variables['hsig'].standard_name='sea_surface_wave_significant_height'
nco.variables['uwind'].standard_name='eastward_wind'
nco.variables['vwind'].standard_name='northward_wind'
# write data
nt=len(jdwave)
uo[:nt,:]=uwind
vo[:nt,:]=vwind
hso[:nt,:]=hs
lato[:]=lati
lono[:]=loni
timeo[:nt]=netCDF4.date2num(jdwave,units=times.units)
nco.close()
nc2=netCDF4.Dataset(urls['FVCOM Currents'])
time_var = nc2.variables['time']
start = jdwave[0]
start = dt.datetime(2005,1,0,0)
# Get desired time step
itime0 = netCDF4.date2index(start,time_var,select='nearest')
dtime = netCDF4.num2date(time_var[itime0],time_var.units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
print 'start: %s' % daystr
stop = jdwave[-1]
# Get desired time step
itime1 = netCDF4.date2index(stop,time_var,select='nearest')
dtime = netCDF4.num2date(time_var[itime1],time_var.units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
print 'stop : %s' % daystr
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-20-d77d6d90f2e8> in <module>() 1 start = jdwave[0] ----> 2 start = dt.datetime(2005,1,0,0) 3 # Get desired time step 4 itime0 = netCDF4.date2index(start,time_var,select='nearest') 5 dtime = netCDF4.num2date(time_var[itime0],time_var.units) ValueError: day is out of range for month
if 'lonc2' not in locals():
lonc2=nc2.variables['lonc'][:]
latc2=nc2.variables['latc'][:]
ind = nearxy(lonc2,latc2,loni,lati)
ind
ind=array([43629, 50361])
#itime0=1000
#itime1=1200
jdvel = netCDF4.num2date(time_var[itime0:itime1],time_var.units)
print jdvel[0].strftime('%Y-%b-%d %H:%M')
def fvcom_adcp(nc,ind,u='u',v='v',istart=0,istop=1,chunk=100):
# single index of element, returns time series from all depths
nz=shape(nc.variables[u])[1]
nt=len(range(istart,istop))
uc=ones((nt,nz))
vc=ones((nt,nz))
time1=time.time()
for i in range(istart,istop,chunk):
itime=arange(i,min(i+chunk,istop))
ichunk = itime-istart
uc[ichunk,:] = nc.variables[u][itime,:,ind]
vc[ichunk,:] = nc.variables[v][itime,:,ind]
time2=time.time()
print 'elapsed seconds = %.1f, %.1f percent done' % (time2-time1,((float(i-istart)+chunk)/nt*100.0))
return uc,vc
nsta=len(ind)
nz=shape(nc2.variables['u'])[1]
nt=len(range(itime0,itime1))
uc = vc = zeros([nt,nz,nsta])
zc = zeros([nz,nsta])
def fvcom_depth_at_face(nc,iface):
# nc = FVCOM Dataset
# ind = single face index (e.g. for currents). Not node index!
nv=nc.variables['nv'][:,iface]
siglay=nc2.variables['siglay'][:,nv]
h = nc.variables['h'][nv]
zc = mean(h*siglay,axis=1)
return zc
for i in range(len(sta)):
uc[:,:,i],vc[:,:,i] = fvcom_adcp(nc2,ind[i],istart=itime0,istop=itime1)
zc[:,i] = fvcom_depth_at_face(nc,ind[i])
uspd=sqrt(uc*uc+vc*vc)
shape(jdvel)
# make the time series plot, with nicely formatted labels
MyDateFormatter = DateFormatter('%b %d\n %Y')
fig = plt.figure(figsize=(22,6), dpi=80)
ax1 = fig.add_subplot(111)
#jd_est = jd - dt.timedelta(hours=5)
# plot entire hindcast/forecast
ax1.plot(jdvel,uspd[:,:,0]) # convert from meters to feet
# plot only more recent levels (from step 72 on...)
#ax1.plot(jd_est[72:],z[72:,:]/.3048) # convert from meters to feet
ax1.xaxis.set_major_formatter(MyDateFormatter)
ax1.grid(True)
setp(gca().get_xticklabels(), rotation=0, horizontalalignment='center')
#plt.title('NECOFS Forecast Waves from %s Forecast' % model)
ax1.set_ylabel('m/s')
ax1.legend(names,loc='upper left')
print nc2.variables['u']
nco = netCDF4.Dataset('/usgs/data1/notebook/kemp_vel.nc', 'w',clobber=True)
nco.createDimension('depth', nz)
nco.createDimension('time', None)
nco.createDimension('station',nsta)
timeo = nco.createVariable('time', 'f4', ('time'))
lono = nco.createVariable('lon', 'f4', ('station'))
lato = nco.createVariable('lat', 'f4', ('station'))
depo = nco.createVariable('depth', 'f4', ('depth','station'))
uo = nco.createVariable('u', 'f4', ('time', 'depth', 'station'))
vo = nco.createVariable('v', 'f4', ('time', 'depth', 'station'))
nco.variables['time'].units = nc2.variables['time'].units
nco.variables['depth'].units = 'm'
nco.variables['lon'].units=nc2.variables['lonc'].units
nco.variables['lat'].units=nc2.variables['latc'].units
nco.variables['u'].units=nc2.variables['u'].units
nco.variables['v'].units=nc2.variables['v'].units
nco.variables['u'].standard_name=nc2.variables['u'].standard_name
nco.variables['v'].standard_name=nc2.variables['v'].standard_name
nt=len(jdvel)
# write data
uo[:nt,:,:]=uc
vo[:nt,:,:]=vc
lato[:]=lati
lono[:]=loni
depo[:,:]=zc
timeo[:nt]=netCDF4.date2num(jdvel,units=time_var.units)
nco.close()