import cmocean.cm as cm
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import xarray as xr
from salishsea_tools import visualisations as vis
from salishsea_tools import viz_tools
%matplotlib inline
mesh = nc.Dataset('/home/sallen/MEOPAR/grid/mesh_mask202108.nc')
bathy = nc.Dataset('/home/sallen/MEOPAR/grid/bathymetry_202108.nc')
tmask = mesh['tmask']
deptht = mesh['gdept_1d'][0]
year = '2007'
olddir = ''
newdir = 'recon'
def surface_plots(tracers, file, month, cmap, cdiff, vmax, vmin, dvmax, twoD=False, olddir=olddir, d201905=False, zoom=[0, 0, 0, 0], dl=0, dosum=False):
fig, axs = plt.subplots(1, 3, figsize=(15, 7))
if d201905:
dm = xr.open_dataset(f'/results/SalishSea/month-avg.201905/SalishSea_1m_{year}{month}_{year}{month}_{file}_T.nc')
else:
dm = xr.open_dataset(
f'/data/sallen/results/MEOPAR/202111/{olddir}/SalishSea_1h_20070101_20070131_{file}_T_20070115-20070115.nc')
if dosum:
dm['summed'] = 0.*dm[tracers[0]]
for tr in tracers:
dm['summed'] = dm['summed'] + dm[tr]
tracer = 'summed'
else:
tracer = tracers
if twoD:
mean_field = np.ma.array(dm[tracer][0], mask=1-tmask[0, 0])
else:
mean_field = np.ma.array(dm[tracer][0, dl], mask=1-tmask[0, dl])
colours = axs[0].pcolormesh(mean_field, cmap=cmap, vmax=vmax, vmin=vmin)
axs[0].set_title(olddir)
fig.colorbar(colours, ax=axs[0])
ds = xr.open_dataset(
f'/data/sallen/results/MEOPAR/202111/{newdir}/SalishSea_1h_20070101_20070115_{file}_T_20070115-20070115.nc')
if dosum:
ds['summed'] = 0.*ds[tracers[0]]
for tr in tracers:
ds['summed'] = ds['summed'] + ds[tr]
if twoD:
field = np.ma.array(ds[tracer][0], mask=1-tmask[0,0])
else:
field = np.ma.array(ds[tracer][0, dl], mask=1-tmask[0, dl])
colours = axs[1].pcolormesh(field, cmap=cmap, vmax=vmax, vmin=vmin)
axs[1].set_title(newdir)
fig.colorbar(colours, ax=axs[1])
colours = axs[2].pcolormesh(field-mean_field, cmap=cdiff, vmax=dvmax, vmin=-dvmax)
axs[2].set_title("New - Old")
fig.colorbar(colours, ax=axs[2])
ds.close()
dm.close()
for ax in axs:
viz_tools.set_aspect(ax)
if sum(zoom) > 0:
ax.set_ylim(zoom[0], zoom[1])
ax.set_xlim(zoom[2], zoom[3])
return fig
def profiles(tracers, file, olddir=olddir, d201905=False, dosum=False):
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
if d201905:
dm = xr.open_dataset(f'/results/SalishSea/month-avg.201905/SalishSea_1m_{year}{month}_{year}{month}_{file}_T.nc')
else:
dm = xr.open_dataset(
f'/data/sallen/results/MEOPAR/202111/{olddir}/SalishSea_1h_20070101_20070131_{file}_T_20070115-20070115.nc')
if dosum:
dm['summed'] = 0.*dm[tracers[0]]
for tr in tracers:
dm['summed'] = dm['summed'] + dm[tr]
tracer = 'summed'
else:
tracer = tracers
mean_field = np.ma.array(dm[tracer][0], mask=1-tmask[0])
axs[0].plot(mean_field.mean(axis=1).mean(axis=1), deptht, linewidth=2, label='Old')
axs[1].plot(np.zeros_like(deptht), deptht, linewidth=2, label='Old')
ds = xr.open_dataset(
f'/data/sallen/results/MEOPAR/202111/{newdir}/SalishSea_1h_20070101_20070115_{file}_T_20070115-20070115.nc')
if dosum:
ds['summed'] = 0.*ds[tracers[0]]
for tr in tracers:
ds['summed'] = ds['summed'] + ds[tr]
field = np.ma.array(ds[tracer][0], mask=1-tmask[0])
axs[0].plot(field.mean(axis=1).mean(axis=1), deptht, label='New')
axs[1].plot((field-mean_field).mean(axis=1).mean(axis=1), deptht, label='New')
print (mean_field.mean(axis=1).mean(axis=1).max(), field.mean(axis=1).mean(axis=1).max())
ds.close()
dm.close()
axs[0].invert_yaxis()
axs[0].legend(loc='best')
axs[1].invert_yaxis()
axs[1].legend(loc='best')
return fig
def thalweg_plots(tracers, file, cmap, cdiff, vmax, vmin, dvmax, olddir=olddir, d201905=False, maxdepth=0, dosum=False):
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
if d201905:
dm = xr.open_dataset(f'/results/SalishSea/month-avg.201905/SalishSea_1m_{year}{month}_{year}{month}_{file}_T.nc')
else:
dm = xr.open_dataset(
f'/data/sallen/results/MEOPAR/202111/{olddir}/SalishSea_1h_20070101_20070131_{file}_T_20070115-20070115.nc')
if dosum:
dm['summed'] = 0.*dm[tracers[0]]
for tr in tracers:
dm['summed'] = dm['summed'] + dm[tr]
tracer = 'summed'
else:
tracer = tracers
mean_field = np.array(dm[tracer][0])
colours = vis.contour_thalweg(axs[0], mean_field, bathy, mesh, np.arange(vmin, 1.1*vmax+0.1*vmin, (vmax-vmin)/10.), cmap=cmap)
axs[0].set_title(olddir)
ds = xr.open_dataset(
f'/data/sallen/results/MEOPAR/202111/{newdir}/SalishSea_1h_20070101_20070115_{file}_T_20070115-20070115.nc')
if dosum:
ds['summed'] = 0.*ds[tracers[0]]
for tr in tracers:
ds['summed'] = ds['summed'] + ds[tr]
field = np.array(ds[tracer][0])
colours = vis.contour_thalweg(axs[1], field, bathy, mesh,np.arange(vmin, 1.1*vmax+0.1*vmin, (vmax-vmin)/10.), cmap=cmap)
axs[1].set_title(newdir)
colours = vis.contour_thalweg(axs[2], field-mean_field, bathy, mesh, np.arange(-dvmax, 1.2*dvmax, dvmax/5),
cmap=cdiff)
axs[2].set_title('New - Old')
if maxdepth > 0:
for ax in axs:
ax.set_ylim(maxdepth, 0)
ds.close()
dm.close()
fig.tight_layout()
return fig
tracer = 'votemper'
file = 'grid'
month = '01'
cmap = cm.thermal
cdiff = cm.balance
vmax = 10
vmin = 4
dvmax = 2
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
9.23042912824859 9.230887276813233
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax, maxdepth=15)
tracer = 'nitrate'
file = 'ptrc'
month = '02'
cmap = cm.rain
cdiff = cm.balance
vmax = 30
vmin = 5
dvmax = 5
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
30.117854110346975 30.119179905610093
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'diatoms'
file = 'ptrc'
month = '03'
cmap = cm.algae
cdiff = cm.balance
vmax = 1
vmin = 0
dvmax = 0.2
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
0.11132826308211242 0.14076791415773973
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'microzooplankton'
file = 'ptrc'
month = '04'
cmap = cm.matter
cmap.set_bad('#8b7765')
cb1 = "Microzooplankton (uM N)"
cdiff = cm.balance
cdiff.set_bad('#8b7765')
cb2 = 'Microzooplankton (uM N)'
vmax = 1
vmin = 0
dvmax = 0.2
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
0.19257068767311925 0.24216836822571358
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'vosaline'
file = 'grid'
month = '05'
cmap = cm.haline
cdiff = cm.balance
vmax = 32
vmin = 15
dvmax = 5
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
31.879404558304984 31.87795229863348
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'flagellates'
file = 'ptrc'
month = '06'
cmap = cm.algae
cdiff = cm.balance
vmax = 1
vmin = 0
dvmax = 0.5
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
0.22420064332130854 0.40474390260024073
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'ammonium'
file = 'ptrc'
month = '07'
cmap = cm.speed
cdiff = cm.balance
vmax = 1
vmin = 0
dvmax = 0.2
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
1.041049437208646 1.0177258493041514
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'sossheig'
file = 'grid'
month = '08'
cmap = cm.tarn
cdiff = cm.balance
vmax = 2
vmin = -2
dvmax = 0.02
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax, twoD=True)
tracer = 'particulate_organic_nitrogen'
file = 'ptrc'
month = '09'
cmap = cm.matter
cdiff = cm.balance
vmax = 0.5
vmin = 0
dvmax = 0.05
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
0.023920590529586328 0.04272459476886933
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'silicon'
file = 'ptrc'
month = '10'
cmap = cm.turbid
cdiff = cm.balance
vmax = 60
vmin = 20
dvmax = 10
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
70.7864824239477 70.7417580509131
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'dissolved_inorganic_carbon'
file = 'carp'
month = '11'
cmap = 'cividis_r'
cdiff = cm.balance
vmax = 2200
vmin = 1800
dvmax = 100
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
2158.6171481949586 2158.5659758901365
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'total_alkalinity'
file = 'carp'
month = '12'
cmap = cm.ice_r
cmap.set_bad('#8b7765')
cdiff = cm.balance
vmax = 2300
vmin = 1900
dvmax = 100
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
2191.077701460487 2190.9661372789205
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
tracer = 'Fraser_tracer'
file = 'carp'
month = '12'
cmap = cm.amp
cdiff = cm.balance
vmax = 20
vmin = 0
dvmax = 2
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax, dl=3,
zoom=[320, 500, 290, 398])
tracer = 'dissolved_oxygen'
file = 'carp'
month = '12'
cmap = cm.oxy
cmap.set_bad('#8b7765')
cdiff = cm.balance
vmax = 250
vmin = 50
dvmax = 50
fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax)
fig = profiles(tracer, file)
258.6935390008605 259.5824079658322
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)
o_gemlam = xr.open_dataset('/results/forcing/atmospheric/GEM2.5/gemlam/gemlam_y2007m01d14.nc')
r_gemlam = xr.open_dataset('/ocean/arandhawa/reconstructed_data_2007/recon_y2007m01d14.nc')
fig, axs = plt.subplots(1, 3, figsize=(20, 5))
o_gemlam.solar.mean(axis=0).plot(ax=axs[0])
r_gemlam.solar.mean(axis=0).plot(ax=axs[1])
(r_gemlam.solar.mean(axis=0)-o_gemlam.solar.mean(axis=0)).plot(ax=axs[2]);
fig, axs = plt.subplots(1, 3, figsize=(20, 5))
o_gemlam.therm_rad.mean(axis=0).plot(ax=axs[0])
r_gemlam.therm_rad.mean(axis=0).plot(ax=axs[1])
(r_gemlam.therm_rad.mean(axis=0)-o_gemlam.therm_rad.mean(axis=0)).plot(ax=axs[2]);
fig, axs = plt.subplots(1, 3, figsize=(20, 5))
o_gemlam.qair.mean(axis=0).plot(ax=axs[0])
r_gemlam.qair.mean(axis=0).plot(ax=axs[1])
(r_gemlam.qair.mean(axis=0)-o_gemlam.qair.mean(axis=0)).plot(ax=axs[2]);
fig, axs = plt.subplots(1, 3, figsize=(20, 5))
o_gemlam.tair.mean(axis=0).plot(ax=axs[0])
r_gemlam.tair.mean(axis=0).plot(ax=axs[1])
(r_gemlam.tair.mean(axis=0)-o_gemlam.tair.mean(axis=0)).plot(ax=axs[2]);
fig, axs = plt.subplots(1, 3, figsize=(20, 5))
o_gemlam.precip.mean(axis=0).plot(ax=axs[0])
r_gemlam.precip.mean(axis=0).plot(ax=axs[1])
(r_gemlam.precip.mean(axis=0)-o_gemlam.precip.mean(axis=0)).plot(ax=axs[2]);
fig, axs = plt.subplots(1, 3, figsize=(20, 5))
o_gemlam.u_wind.mean(axis=0).plot(ax=axs[0])
r_gemlam.u_wind.mean(axis=0).plot(ax=axs[1])
(r_gemlam.u_wind.mean(axis=0)-o_gemlam.u_wind.mean(axis=0)).plot(ax=axs[2]);
fig, axs = plt.subplots(1, 3, figsize=(20, 5))
o_gemlam.v_wind.mean(axis=0).plot(ax=axs[0])
r_gemlam.v_wind.mean(axis=0).plot(ax=axs[1])
(r_gemlam.v_wind.mean(axis=0)-o_gemlam.v_wind.mean(axis=0)).plot(ax=axs[2]);