In [1]:
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from salishsea_tools import geo_tools, viz_tools, evaltools as et
import gsw
import netCDF4 as nc
import datetime as dt
import cmocean
import pickle
import os
import warnings
warnings.filterwarnings("ignore")
%matplotlib inline
In [2]:
obslat=50.117
obslon=-125.223
In [3]:
with nc.Dataset('/ocean/eolson/MEOPAR/NEMO-forcing/grid/mesh_mask201702_noLPE.nc') as mesh:
tmask=np.copy(mesh.variables['tmask'][0,:,:,:])
navlat=np.copy(mesh.variables['nav_lat'][:,:])
navlon=np.copy(mesh.variables['nav_lon'][:,:])
e3t_0=np.copy(mesh.variables['e3t_0'][0,:,:,:])
print(mesh.variables.keys())
dict_keys(['nav_lon', 'nav_lat', 'nav_lev', 'time_counter', 'tmask', 'umask', 'vmask', 'fmask', 'tmaskutil', 'umaskutil', 'vmaskutil', 'fmaskutil', 'glamt', 'glamu', 'glamv', 'glamf', 'gphit', 'gphiu', 'gphiv', 'gphif', 'e1t', 'e1u', 'e1v', 'e1f', 'e2t', 'e2u', 'e2v', 'e2f', 'ff', 'mbathy', 'misf', 'isfdraft', 'e3t_0', 'e3u_0', 'e3v_0', 'e3w_0', 'gdept_0', 'gdepu', 'gdepv', 'gdepw_0', 'gdept_1d', 'gdepw_1d', 'e3t_1d', 'e3w_1d'])
In [4]:
j,i=geo_tools.find_closest_model_point(obslon,obslat,navlon,navlat)
fig,ax=plt.subplots(1,2,figsize=(10,5))
ax[0].contour(navlon,navlat,tmask[0,:,:],(0.5,))
ax[0].plot(navlon[j,i],navlat[j,i],'r*')
ax[0].set_xlim(-125.4,-125)
ax[0].set_ylim(50,50.4)
viz_tools.set_aspect(ax[0],coords='map')
ax[0].set_title('Geo coords')
ax[1].contour(tmask[0,:,:],(.5,))
ax[1].plot(i,j,'r*')
ax[1].set_xlim(100,180)
ax[1].set_ylim(730,800)
viz_tools.set_aspect(ax[1],coords='grid')
ax[1].set_title('Grid coords')
ks=np.sum(tmask[:,j,i])
print('depth is {0} grid points at site or roughly {1:3.1f} m'.format(ks,np.sum(e3t_0[:ks,j,i])))
depth is 15 grid points at site or roughly 15.1 m
In [5]:
print(j,i)
761 137
In [5]:
df=pd.read_csv('/ocean/eolson/MEOPAR/obs/IzettCTD/Quadra_Hakai_BoL-data_201906-202004.txt',header=0,skiprows=(1,))
In [6]:
df.head()
Out[6]:
| Days | TSG_T | TSG_S | CTD_1_T | CTD_1_S | CTD_1_pr | CTD_2_T | CTD_2_S | CTD_2_pr | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 352.144785 | 7.259759 | 21.940098 | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 352.148303 | 7.269010 | 21.991546 | NaN | NaN | NaN | NaN | NaN | NaN |
| 2 | 352.151799 | 7.278374 | 22.002477 | NaN | NaN | NaN | NaN | NaN | NaN |
| 3 | 352.155294 | 7.269549 | 21.972746 | NaN | NaN | NaN | NaN | NaN | NaN |
| 4 | 352.158790 | 7.282700 | 22.012073 | NaN | NaN | NaN | NaN | NaN | NaN |
In [7]:
df['dtUTC']=[dt.datetime(2014,1,1)+dt.timedelta(days=ii) for ii in df['Days']]
In [8]:
df['SA']=[gsw.SA_from_SP(isp,1.4,obslon,obslat) for isp in df['TSG_S']]
In [9]:
df['CT']=[gsw.CT_from_t(isa,ist,1.4) for isa,ist in zip(df['SA'],df['TSG_T'])]
In [10]:
df['Z']=1.4
df['Lat']=obslat
df['Lon']=obslon
In [11]:
df['k']=1
df['i']=i
df['j']=j
In [12]:
df.head()
Out[12]:
| Days | TSG_T | TSG_S | CTD_1_T | CTD_1_S | CTD_1_pr | CTD_2_T | CTD_2_S | CTD_2_pr | dtUTC | SA | CT | Z | Lat | Lon | k | i | j | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 352.144785 | 7.259759 | 21.940098 | NaN | NaN | NaN | NaN | NaN | NaN | 2014-12-19 03:28:29.424000 | 22.046300 | 7.429968 | 1.4 | 50.117 | -125.223 | 1 | 137 | 761 |
| 1 | 352.148303 | 7.269010 | 21.991546 | NaN | NaN | NaN | NaN | NaN | NaN | 2014-12-19 03:33:33.379200 | 22.097997 | 7.438737 | 1.4 | 50.117 | -125.223 | 1 | 137 | 761 |
| 2 | 352.151799 | 7.278374 | 22.002477 | NaN | NaN | NaN | NaN | NaN | NaN | 2014-12-19 03:38:35.433600 | 22.108981 | 7.448129 | 1.4 | 50.117 | -125.223 | 1 | 137 | 761 |
| 3 | 352.155294 | 7.269549 | 21.972746 | NaN | NaN | NaN | NaN | NaN | NaN | 2014-12-19 03:43:37.401600 | 22.079106 | 7.439522 | 1.4 | 50.117 | -125.223 | 1 | 137 | 761 |
| 4 | 352.158790 | 7.282700 | 22.012073 | NaN | NaN | NaN | NaN | NaN | NaN | 2014-12-19 03:48:39.456000 | 22.118623 | 7.452410 | 1.4 | 50.117 | -125.223 | 1 | 137 | 761 |
In [13]:
flen=1
namfmt='nowcast'
filemap={'vosaline':'grid_T','votemper':'grid_T'}
fdict={'grid_T':1,}
PATH='/results2/SalishSea/nowcast-green.201905/'
runID='HC201905'
dstart=dt.datetime(2014,12,19)
dend=dt.datetime(2020,6,20)
In [14]:
len(df)
Out[14]:
450137
In [15]:
# remove a day because it is currently missing from model output
df.drop(df.loc[(df.dtUTC>=dt.datetime(2020,6,10))&(df.dtUTC<dt.datetime(2020,6,11))].index,inplace=True)
In [16]:
data=et.matchData(df,filemap,fdict,dstart,dend,namfmt,PATH,flen,preIndexed=True)
progress: 0.0% progress: 1.1115063133558598% progress: 2.2230126267117196% progress: 3.3345189400675794% progress: 4.446025253423439% progress: 5.557531566779299% progress: 6.669037880135159% progress: 7.7805441934910196% progress: 8.892050506846878% progress: 10.00355682020274% progress: 11.115063133558598% progress: 12.226569446914457% progress: 13.338075760270318% progress: 14.449582073626177% progress: 15.561088386982039% progress: 16.672594700337896% progress: 17.784101013693757% progress: 18.895607327049618% progress: 20.00711364040548% progress: 21.118619953761335% progress: 22.230126267117196% progress: 23.341632580473057% progress: 24.453138893828914% progress: 25.56464520718478% progress: 26.676151520540635% progress: 27.787657833896496% progress: 28.899164147252353% progress: 30.010670460608218% progress: 31.122176773964078% progress: 32.23368308731994% progress: 33.34518940067579% progress: 34.45669571403165% progress: 35.568202027387514% progress: 36.679708340743375% progress: 37.791214654099235% progress: 38.902720967455096% progress: 40.01422728081096% progress: 41.12573359416681% progress: 42.23723990752267% progress: 43.34874622087854% progress: 44.46025253423439% progress: 45.57175884759025% progress: 46.683265160946114% progress: 47.794771474301974% progress: 48.90627778765783% progress: 50.017784101013696% progress: 51.12929041436956% progress: 52.24079672772541% progress: 53.35230304108127% progress: 54.46380935443713% progress: 55.57531566779299% progress: 56.68682198114885% progress: 57.79832829450471% progress: 58.90983460786057% progress: 60.021340921216435% progress: 61.132847234572296% progress: 62.244353547928156% progress: 63.35585986128402% progress: 64.46736617463988% progress: 65.57887248799574% progress: 66.69037880135159% progress: 67.80188511470745% progress: 68.9133914280633% progress: 70.02489774141917% progress: 71.13640405477503% progress: 72.2479103681309% progress: 73.35941668148675% progress: 74.47092299484261% progress: 75.58242930819847% progress: 76.69393562155433% progress: 77.80544193491019% progress: 78.91694824826605% progress: 80.02845456162191% progress: 81.13996087497777% progress: 82.25146718833362% progress: 83.36297350168948% progress: 84.47447981504534% progress: 85.58598612840122% progress: 86.69749244175708% progress: 87.80899875511294% progress: 88.92050506846878% progress: 90.03201138182465% progress: 91.1435176951805% progress: 92.25502400853637% progress: 93.36653032189223% progress: 94.47803663524809% progress: 95.58954294860395% progress: 96.70104926195981% progress: 97.81255557531566% progress: 98.92406188867153%
In [17]:
fig,ax=plt.subplots(1,2,figsize=(12,4))
_,_,_,m0=ax[0].hist2d(data['SA'],data['mod_vosaline'],bins=np.linspace(21,28,100),
cmap=cmocean.cm.amp,cmin=1,cmax=800,vmin=1,vmax=350)
fig.colorbar(m0,ax=ax[0],extend='max')
ax[0].set_aspect(1)
ax[0].set_xlim(21,28)
ax[0].set_ylim(21,28)
ax[0].plot((21,28),(21,28),'k-',lw=.5);
ax[0].set_xlabel('Obs Abs. Sal (g kg$^{-1}$)')
ax[0].set_ylabel('Model Abs. Sal (g kg$^{-1}$)')
_,_,_,m1=ax[1].hist2d(data['CT'],data['mod_votemper'],bins=np.linspace(5.5,10,100),
cmap=cmocean.cm.amp,cmin=1,cmax=350)
fig.colorbar(m1,ax=ax[1])
ax[1].set_aspect(1)
ax[1].set_xlim(5.5,10)
ax[1].set_ylim(5.5,10)
ax[1].plot((5.5,10),(5.5,10),'k-',lw=.5)
ax[1].set_xlabel('Obs Cons. Temp. (deg C)')
ax[1].set_ylabel('Model Cons. Temp. (deg C)')
Out[17]:
Text(0, 0.5, 'Model Cons. Temp. (deg C)')
In [18]:
s_sa=et.stats(data['SA'],data['mod_vosaline'])
print('Absolute Salinity\n N: {0}\n bias: {3}\n RMSE: {4}\n Willmott Skill Score: {5}'.format(*s_sa))
Absolute Salinity N: 447830 bias: 0.8393380791708118 RMSE: 1.4592866685270445 Willmott Skill Score: 0.7603133483217335
In [19]:
s_ct=et.stats(data['CT'],data['mod_votemper'])
print('Conservative Temperature\n N: {0}\n bias: {3}\n RMSE: {4}\n Willmott Skill Score: {5}'.format(*s_ct))
Conservative Temperature N: 447830 bias: -0.8051871415688563 RMSE: 1.2563920299941473 Willmott Skill Score: 0.9780175111340039
In [20]:
df2=data.copy(deep=True)
df2['k']=0
In [21]:
data2=et.matchData(df2,filemap,fdict,dstart,dend,namfmt,PATH,flen,preIndexed=True)
progress: 0.0% progress: 1.1115063133558598% progress: 2.2230126267117196% progress: 3.3345189400675794% progress: 4.446025253423439% progress: 5.557531566779299% progress: 6.669037880135159% progress: 7.7805441934910196% progress: 8.892050506846878% progress: 10.00355682020274% progress: 11.115063133558598% progress: 12.226569446914457% progress: 13.338075760270318% progress: 14.449582073626177% progress: 15.561088386982039% progress: 16.672594700337896% progress: 17.784101013693757% progress: 18.895607327049618% progress: 20.00711364040548% progress: 21.118619953761335% progress: 22.230126267117196% progress: 23.341632580473057% progress: 24.453138893828914% progress: 25.56464520718478% progress: 26.676151520540635% progress: 27.787657833896496% progress: 28.899164147252353% progress: 30.010670460608218% progress: 31.122176773964078% progress: 32.23368308731994% progress: 33.34518940067579% progress: 34.45669571403165% progress: 35.568202027387514% progress: 36.679708340743375% progress: 37.791214654099235% progress: 38.902720967455096% progress: 40.01422728081096% progress: 41.12573359416681% progress: 42.23723990752267% progress: 43.34874622087854% progress: 44.46025253423439% progress: 45.57175884759025% progress: 46.683265160946114% progress: 47.794771474301974% progress: 48.90627778765783% progress: 50.017784101013696% progress: 51.12929041436956% progress: 52.24079672772541% progress: 53.35230304108127% progress: 54.46380935443713% progress: 55.57531566779299% progress: 56.68682198114885% progress: 57.79832829450471% progress: 58.90983460786057% progress: 60.021340921216435% progress: 61.132847234572296% progress: 62.244353547928156% progress: 63.35585986128402% progress: 64.46736617463988% progress: 65.57887248799574% progress: 66.69037880135159% progress: 67.80188511470745% progress: 68.9133914280633% progress: 70.02489774141917% progress: 71.13640405477503% progress: 72.2479103681309% progress: 73.35941668148675% progress: 74.47092299484261% progress: 75.58242930819847% progress: 76.69393562155433% progress: 77.80544193491019% progress: 78.91694824826605% progress: 80.02845456162191% progress: 81.13996087497777% progress: 82.25146718833362% progress: 83.36297350168948% progress: 84.47447981504534% progress: 85.58598612840122% progress: 86.69749244175708% progress: 87.80899875511294% progress: 88.92050506846878% progress: 90.03201138182465% progress: 91.1435176951805% progress: 92.25502400853637% progress: 93.36653032189223% progress: 94.47803663524809% progress: 95.58954294860395% progress: 96.70104926195981% progress: 97.81255557531566% progress: 98.92406188867153%
In [22]:
fig,ax=plt.subplots(1,2,figsize=(12,4))
_,_,_,m0=ax[0].hist2d(data2['SA'],data2['mod_vosaline'],bins=np.linspace(21,28,100),
cmap=cmocean.cm.amp,cmin=1,cmax=800,vmin=1,vmax=350)
fig.colorbar(m0,ax=ax[0],extend='max')
ax[0].set_aspect(1)
ax[0].set_xlim(21,28)
ax[0].set_ylim(21,28)
ax[0].plot((21,28),(21,28),'k-',lw=.5);
ax[0].set_xlabel('Obs Abs. Sal (g kg$^{-1}$)')
ax[0].set_ylabel('Model Abs. Sal (g kg$^{-1}$)')
_,_,_,m1=ax[1].hist2d(data2['CT'],data2['mod_votemper'],bins=np.linspace(5.5,10,100),
cmap=cmocean.cm.amp,cmin=1,cmax=350)
fig.colorbar(m1,ax=ax[1])
ax[1].set_aspect(1)
ax[1].set_xlim(5.5,10)
ax[1].set_ylim(5.5,10)
ax[1].plot((5.5,10),(5.5,10),'k-',lw=.5)
ax[1].set_xlabel('Obs Cons. Temp. (deg C)')
ax[1].set_ylabel('Model Cons. Temp. (deg C)')
Out[22]:
Text(0, 0.5, 'Model Cons. Temp. (deg C)')
In [23]:
s_sa=et.stats(data2['SA'],data2['mod_vosaline'])
print('Absolute Salinity\n N: {0}\n bias: {3}\n RMSE: {4}\n Willmott Skill Score: {5}'.format(*s_sa))
Absolute Salinity N: 447830 bias: 0.7309802323394052 RMSE: 1.3915637204523503 Willmott Skill Score: 0.7874123043385386
In [24]:
s_ct=et.stats(data2['CT'],data2['mod_votemper'])
print('Conservative Temperature\n N: {0}\n bias: {3}\n RMSE: {4}\n Willmott Skill Score: {5}'.format(*s_ct))
Conservative Temperature N: 447830 bias: -0.5625343779442051 RMSE: 1.0545691219621802 Willmott Skill Score: 0.9854292915737839
In [25]:
fig,ax=plt.subplots(2,1,figsize=(16,4))
fig.subplots_adjust(hspace=.3)
ax[0].plot(data['dtUTC'],data['mod_votemper'],'r.',label='model',ms=1)
ax[0].plot(data['dtUTC'],data['CT'],'k.',label='obs',ms=1)
ax[0].legend()
ax[0].set_ylabel('Cons. Temp.\n (deg C)')
ax[1].plot(data['dtUTC'],data['mod_vosaline'],'r.',label='model',ms=1)
ax[1].plot(data['dtUTC'],data['SA'],'k.',label='obs',ms=1)
ax[1].legend()
ax[1].set_ylabel('Abs. Sal.\n (g kg$^{-1}$)')
Out[25]:
Text(0, 0.5, 'Abs. Sal.\n (g kg$^{-1}$)')
SSH Comparisons¶
In [ ]:
In [26]:
dfssh=pd.read_csv('/ocean/eolson/MEOPAR/obs/IzettCTD/WhaletownTidePred/8038-01-JAN-2015_slev.csv',
skiprows=7,delimiter=',',index_col=None,usecols=[0,1])
In [27]:
dfssh
Out[27]:
| Obs_date | SLEV(metres) | |
|---|---|---|
| 0 | 2017/06/16 22:00 | 2.46 |
| 1 | 2017/06/16 23:00 | 2.19 |
| 2 | 2017/06/17 00:00 | 2.11 |
| 3 | 2017/06/17 01:00 | 2.23 |
| 4 | 2017/06/17 02:00 | 2.61 |
| ... | ... | ... |
| 13652 | 2019/06/13 20:00 | 2.79 |
| 13653 | 2019/06/13 21:00 | 3.42 |
| 13654 | 2019/06/13 22:00 | 3.91 |
| 13655 | 2019/06/13 23:00 | 4.15 |
| 13656 | 2019/06/14 00:00 | 4.17 |
13657 rows × 2 columns
In [28]:
dfssh['dtUTC']=[dt.datetime.strptime(ii, '%Y/%m/%d %H:%M') for ii in dfssh['Obs_date']]
In [29]:
jssh,issh=geo_tools.find_closest_model_point(-125.052,50.108,navlon,navlat)
# point is on land so move by 1
issh=issh-1
dfssh['i']=issh
dfssh['j']=jssh
In [30]:
dfssh
Out[30]:
| Obs_date | SLEV(metres) | dtUTC | i | j | |
|---|---|---|---|---|---|
| 0 | 2017/06/16 22:00 | 2.46 | 2017-06-16 22:00:00 | 160 | 747 |
| 1 | 2017/06/16 23:00 | 2.19 | 2017-06-16 23:00:00 | 160 | 747 |
| 2 | 2017/06/17 00:00 | 2.11 | 2017-06-17 00:00:00 | 160 | 747 |
| 3 | 2017/06/17 01:00 | 2.23 | 2017-06-17 01:00:00 | 160 | 747 |
| 4 | 2017/06/17 02:00 | 2.61 | 2017-06-17 02:00:00 | 160 | 747 |
| ... | ... | ... | ... | ... | ... |
| 13652 | 2019/06/13 20:00 | 2.79 | 2019-06-13 20:00:00 | 160 | 747 |
| 13653 | 2019/06/13 21:00 | 3.42 | 2019-06-13 21:00:00 | 160 | 747 |
| 13654 | 2019/06/13 22:00 | 3.91 | 2019-06-13 22:00:00 | 160 | 747 |
| 13655 | 2019/06/13 23:00 | 4.15 | 2019-06-13 23:00:00 | 160 | 747 |
| 13656 | 2019/06/14 00:00 | 4.17 | 2019-06-14 00:00:00 | 160 | 747 |
13657 rows × 5 columns
In [31]:
flen=1
namfmt='nowcast'
filemap={'sossheig':'grid_T',}
fdict={'grid_T':1,}
PATH='/results2/SalishSea/nowcast-green.201905/'
runID='HC201905'
dstart=dt.datetime(2017,6,1)
dend=dt.datetime(2019,6,15)
In [32]:
dataSSH=et.matchData(dfssh,filemap,fdict,dstart,dend,namfmt,PATH,flen,preIndexed=True,sdim=2)
progress: 0.0% progress: 36.61126162407557% progress: 73.22252324815113%
In [33]:
fig,ax=plt.subplots(1,1,figsize=(16,4))
ax.plot(dataSSH['dtUTC'],dataSSH['mod_sossheig']-np.mean(dataSSH['mod_sossheig']),'r.',label='model',ms=1)
ax.plot(dataSSH['dtUTC'],dataSSH['SLEV(metres)']-np.mean(dataSSH['SLEV(metres)']),'k.',label='obs',ms=1)
ax.legend()
ax.set_ylabel('SSH (m)')
Out[33]:
Text(0, 0.5, 'SSH (m)')
In [34]:
fig,ax=plt.subplots(1,1,figsize=(16,4))
ax.plot(dataSSH['dtUTC']+dt.timedelta(hours=.5),dataSSH['mod_sossheig']-np.mean(dataSSH['mod_sossheig']),'r-',label='model',ms=1)
# model is actually value on half hour but was matched to obs on the hour, so added 30 min to model times
ax.plot(dataSSH['dtUTC'],dataSSH['SLEV(metres)']-np.mean(dataSSH['SLEV(metres)']),'k-',label='obs',ms=1)
ax.legend()
ax.set_ylabel('SSH (m)')
ax.set_xlim((dt.datetime(2017,7,1),dt.datetime(2017,7,30)))
Out[34]:
(736511.0, 736540.0)
In [35]:
fig,ax=plt.subplots(1,1,figsize=(16,4))
ax.plot(dataSSH['dtUTC']+dt.timedelta(hours=.5),dataSSH['mod_sossheig']-np.mean(dataSSH['mod_sossheig']),'r-',label='model',ms=1)
ax.plot(dataSSH['dtUTC'],dataSSH['SLEV(metres)']-np.mean(dataSSH['SLEV(metres)']),'k-',label='obs',ms=1)
ax.legend()
ax.set_ylabel('SSH (m)')
ax.set_xlim((dt.datetime(2018,9,1),dt.datetime(2018,9,30)))
Out[35]:
(736938.0, 736967.0)
In [ ]: