Looking at the 2016 renewal in nowcast-green, nowcast and obs
In [1]:
import matplotlib.pyplot as plt
import numpy as np
import netCDF4 as nc
import datetime
from dateutil import tz
import os
import pandas as pd
import xarray as xr
from salishsea_tools import (
geo_tools,
places,
psu_tools,
teos_tools,
data_tools,
tidetools,
)
from nowcast import analyze
from nowcast.figures import shared
%matplotlib inline
In [2]:
import seaborn as sns
sns.set_context('talk')
sns.set_style('darkgrid')
sns.set_color_codes()
In [3]:
runs={}
t_o=datetime.datetime(2016,1,1); t_f = datetime.datetime(2016,10,10)
fnames = analyze.get_filenames(t_o, t_f, '1d', 'grid_T', '/results/SalishSea/nowcast-green/')
grid_B=nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_downonegrid2.nc')
mesh_mask = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/mesh_mask_downbyone2.nc')
runs = {'nowcast-green': {'grid': grid_B,
'mesh': mesh_mask,
'fnames': fnames,
'nemo36': True}}
In [11]:
fnames = analyze.get_filenames(t_o, t_f, '1d', 'grid_T', '/results/SalishSea/nowcast/')
grid_B=nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
mesh_mask = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/mesh_mask_SalishSea2.nc')
runs['nowcast']= {'grid': grid_B,
'mesh': mesh_mask,
'fnames': fnames,
'nemo36': False}
In [4]:
def get_onc_TS_time_series(station, t_o, t_f):
"""Grab the ONC temperature and salinity time series for a station between dates
t_o and t_f. Return results as separate temperature and salinty data frames."""
numdays = (t_f-t_o).days
dates = [t_o + datetime.timedelta(days=num)
for num in range(0, numdays+1)]
sal_pd = pd.DataFrame({'time':[],
'data': []})
temp_pd = pd.DataFrame({'time': [],
'data': []})
station_code = places.PLACES[station]['ONC stationCode']
for date in dates:
onc_data = data_tools.get_onc_data(
'scalardata', 'getByStation', os.environ['ONC_USER_TOKEN'],
station=station_code,
deviceCategory='CTD', sensors='salinity,temperature',
dateFrom=data_tools.onc_datetime(date, 'utc'))
try:
ctd_data=data_tools.onc_json_to_dataset(onc_data, teos=False) #keep in PSU!
#quality control
qc_sal = np.array(ctd_data.salinity.qaqcFlag)
qc_temp = np.array(ctd_data.temperature.qaqcFlag)
#average
sal_pd = sal_pd.append({'time': ctd_data.salinity.sampleTime.values[0],
'data': ctd_data.salinity.values[qc_sal==1].mean()},
ignore_index=True)
temp_pd = temp_pd.append({'time': ctd_data.temperature.sampleTime.values[0],
'data': ctd_data.temperature.values[qc_temp==1].mean()},
ignore_index=True)
except TypeError:
print('No data for {} at {}'.format(date, station))
return sal_pd, temp_pd
In [5]:
def get_model_time_series(station, fnames, grid_B, mesh_mask, nemo_36=True):
"""Retrieve the density, salinity and temperature time series at a station.
Time series is created from files listed in fnames"""
if nemo_36:
depth_var='gdept_0'
depth_var_w = 'gdepw_0'
else:
depth_var='gdept'
depth_var_w = 'gdepw'
#station info
lon = places.PLACES[station]['lon lat'][0]
lat = places.PLACES[station]['lon lat'][1]
depth = places.PLACES[station]['depth']
# model corresponding locations and variables
bathy, X, Y = tidetools.get_bathy_data(grid_B)
j, i = geo_tools.find_closest_model_point(lon,lat,X,Y, land_mask=bathy.mask)
model_depths = mesh_mask.variables[depth_var][0,:,j,i]
tmask = mesh_mask.variables['tmask'][0,:,j,i]
wdeps = mesh_mask.variables[depth_var_w][0,:,j,i]
sal, time = analyze.combine_files(fnames,'vosaline','None',j,i)
temp, time = analyze.combine_files(fnames,'votemper','None',j,i)
# interpolate:
sal_interp=np.array(
[shared.interpolate_tracer_to_depths(
sal[i,:],model_depths,depth,tmask,wdeps)
for i in range(sal.shape[0])])
temp_interp=np.array(
[shared.interpolate_tracer_to_depths(
temp[i,:],model_depths,depth,tmask,wdeps)
for i in range(temp.shape[0])])
# convert to psu for using density function
if nemo_36:
sal_interp = teos_tools.teos_psu(sal_interp)
density = psu_tools.calculate_density(temp_interp, sal_interp)
return density, sal_interp, temp_interp, time
In [6]:
obs_sal={}
obs_temp={}
for station in ['Central node']:
obs_sal[station], obs_temp[station] = get_onc_TS_time_series(station, t_o, t_f)
/home/nsoontie/anaconda3/envs/analysis/lib/python3.5/site-packages/numpy/core/_methods.py:59: RuntimeWarning: Mean of empty slice.
warnings.warn("Mean of empty slice.", RuntimeWarning)
No data for 2016-05-02 00:00:00 at Central node No data for 2016-08-15 00:00:00 at Central node
/home/nsoontie/anaconda3/envs/analysis/lib/python3.5/site-packages/ipykernel/__main__.py:28: VisibleDeprecationWarning: boolean index did not match indexed array along dimension 0; dimension is 86398 but corresponding boolean dimension is 86397
No data for 2016-10-08 00:00:00 at Central node No data for 2016-10-09 00:00:00 at Central node No data for 2016-10-10 00:00:00 at Central node
In [12]:
rhos = {'nowcast': {}, 'nowcast-green': {}}
times = {'nowcast': {}, 'nowcast-green': {}}
sals={'nowcast': {}, 'nowcast-green': {}}
temps={'nowcast': {}, 'nowcast-green': {}}
stations = ['Central node',]
for sim in ['nowcast', 'nowcast-green']:
print(sim)
for station in stations:
print(station)
rhos[sim][station], sals[sim][station], temps[sim][station], times[sim][station] = \
get_model_time_series(station, runs[sim]['fnames'], runs[sim]['grid'],
runs[sim]['mesh'], nemo_36=runs[sim]['nemo36'] )
nowcast Central node nowcast-green Central node
In [10]:
runs.keys()
Out[10]:
dict_keys(['nowcast-green'])
In [22]:
fig, axs = plt.subplots(1,1,figsize=(15,5), sharex=True)
names = [ 'Salinty [g/kg]',]
titles = ['salinity',]
ticks = [[29.6, 32],]
cols = ['r','b']
labels={'nowcast-green': 'Version 2', 'nowcast': 'Version 1'}
#obs=['',]
obs = [obs_sal,]
for i, station in enumerate(stations):
axc = [axs,]
for sim, col in zip(['nowcast-green','nowcast'], cols):
variables = [sals[sim], ]
t = times[sim]
for var, name, title, ax, ob, tick in zip(variables, names, titles, axc, obs,ticks):
if sim == 'nowcast-green': #only plot obs once
ax.plot(ob[station].time, teos_tools.psu_teos(ob[station].data), 'k', label='Observations', lw=2)
ax.plot(np.array(t[station]), teos_tools.psu_teos(np.array(var[station])), c=col,
label=labels[sim], lw=2)
ax.set_ylabel('Daily averaged {}'.format(name))
ax.set_title('{} - {} m'.format(station, places.PLACES[station]['depth']))
ax.set_ylim(tick)
for ax in axc:
ax.get_yaxis().get_major_formatter().set_useOffset(False)
ax.legend(loc=0)
fig.autofmt_xdate()
ax.set_xlim([datetime.datetime(2016,2,14), datetime.datetime(2016,9,25)])
Out[22]:
(736008.0, 736232.0)
In [23]:
fig.savefig('/home/nsoontie/Desktop/Central_2016.png', dpi=300, bbox_inches='tight')
In [ ]: