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 = 'ave3h'
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)
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/file_manager.py in _acquire_with_cache_info(self, needs_lock) 198 try: --> 199 file = self._cache[self._key] 200 except KeyError: /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/lru_cache.py in __getitem__(self, key) 52 with self._lock: ---> 53 value = self._cache[key] 54 self._cache.move_to_end(key) KeyError: [<class 'netCDF4._netCDF4.Dataset'>, ('/data/sallen/results/MEOPAR/202111/base_32/01jan07/SalishSea_1h_20070101_20070131_grid_T_20070115-20070115.nc',), 'r', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False))] During handling of the above exception, another exception occurred: FileNotFoundError Traceback (most recent call last) <ipython-input-12-cfcef9b9753e> in <module> 7 vmin = 4 8 dvmax = 2 ----> 9 fig = surface_plots(tracer, file, month, cmap, cdiff, vmax, vmin, dvmax) <ipython-input-5-08ae36a32687> in surface_plots(tracers, file, month, cmap, cdiff, vmax, vmin, dvmax, twoD, olddir, d201905, zoom, dl, dosum) 5 dm = xr.open_dataset(f'/results/SalishSea/month-avg.201905/SalishSea_1m_{year}{month}_{year}{month}_{file}_T.nc') 6 else: ----> 7 dm = xr.open_dataset( 8 f'/data/sallen/results/MEOPAR/202111/{olddir}/SalishSea_1h_20070101_20070131_{file}_T_20070115-20070115.nc') 9 /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/api.py in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, backend_kwargs, *args, **kwargs) 494 495 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None) --> 496 backend_ds = backend.open_dataset( 497 filename_or_obj, 498 drop_variables=drop_variables, /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/netCDF4_.py in open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, format, clobber, diskless, persist, lock, autoclose) 547 548 filename_or_obj = _normalize_path(filename_or_obj) --> 549 store = NetCDF4DataStore.open( 550 filename_or_obj, 551 mode=mode, /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/netCDF4_.py in open(cls, filename, mode, format, group, clobber, diskless, persist, lock, lock_maker, autoclose) 378 netCDF4.Dataset, filename, mode=mode, kwargs=kwargs 379 ) --> 380 return cls(manager, group=group, mode=mode, lock=lock, autoclose=autoclose) 381 382 def _acquire(self, needs_lock=True): /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/netCDF4_.py in __init__(self, manager, group, mode, lock, autoclose) 326 self._group = group 327 self._mode = mode --> 328 self.format = self.ds.data_model 329 self._filename = self.ds.filepath() 330 self.is_remote = is_remote_uri(self._filename) /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/netCDF4_.py in ds(self) 387 @property 388 def ds(self): --> 389 return self._acquire() 390 391 def open_store_variable(self, name, var): /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/netCDF4_.py in _acquire(self, needs_lock) 381 382 def _acquire(self, needs_lock=True): --> 383 with self._manager.acquire_context(needs_lock) as root: 384 ds = _nc4_require_group(root, self._group, self._mode) 385 return ds /ocean/sallen/miniconda3/envs/py39/lib/python3.9/contextlib.py in __enter__(self) 115 del self.args, self.kwds, self.func 116 try: --> 117 return next(self.gen) 118 except StopIteration: 119 raise RuntimeError("generator didn't yield") from None /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/file_manager.py in acquire_context(self, needs_lock) 185 def acquire_context(self, needs_lock=True): 186 """Context manager for acquiring a file.""" --> 187 file, cached = self._acquire_with_cache_info(needs_lock) 188 try: 189 yield file /ocean/sallen/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/backends/file_manager.py in _acquire_with_cache_info(self, needs_lock) 203 kwargs = kwargs.copy() 204 kwargs["mode"] = self._mode --> 205 file = self._opener(*self._args, **kwargs) 206 if self._mode == "w": 207 # ensure file doesn't get overriden when opened again src/netCDF4/_netCDF4.pyx in netCDF4._netCDF4.Dataset.__init__() src/netCDF4/_netCDF4.pyx in netCDF4._netCDF4._ensure_nc_success() FileNotFoundError: [Errno 2] No such file or directory: b'/data/sallen/results/MEOPAR/202111/base_32/01jan07/SalishSea_1h_20070101_20070131_grid_T_20070115-20070115.nc'
fig = profiles(tracer, file)
9.23042912824859 9.230482687204878
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.117445646507353
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.12316961426497562
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.20687974391729863
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.878966936272686
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.28009593982961734
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.047317910650469
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.028426609683075477
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.7998071505998
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.604579919467
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 2191.0814705201715
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 258.3412151875618
fig = thalweg_plots(tracer, file, cmap, cdiff, vmax, vmin, dvmax)