#!/usr/bin/env python # coding: utf-8 # 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 get_ipython().run_line_magic('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) # 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'] ) # In[10]: runs.keys() # 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)]) # In[23]: fig.savefig('/home/nsoontie/Desktop/Central_2016.png', dpi=300, bbox_inches='tight') # In[ ]: