Notebook

In [2]:

import cmocean.cm as cm
import datetime as dt
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.colors import LogNorm
import numpy as np
import pandas as pd
# import statfieldodels.api as sm
import xarray as xr

from salishsea_tools import viz_tools
import netCDF4 as nc
import arrow
import copy
import math
import matplotlib.dates as mdates

from scipy.io import loadmat

In [3]:

mycmap = copy.copy(cm.oxy)
mycmap.set_bad('darkgreen')

ERDDAP server - https://nwem.apl.washington.edu/erddap/index.html https://nwem.apl.washington.edu/erddap/search/index.html?page=1&itemsPerPage=1000&searchFor=Oxygen - data

In [4]:

with xr.open_dataset('/ocean/atall/MOAD/Obs/PugetSound/ERDDAP/L2_gridded_025_TW_20070101-20071231.nc') as TW_data:
    oxy_TW_obs = TW_data['mass_concentration_of_oxygen_in_sea_water']*44.64
    temp_TW_obs = TW_data['sea_water_temperature']
    salt_TW_obs = TW_data['sea_water_practical_salinity']
    no3_TW_obs = TW_data['mole_concentration_of_nitrate_in_sea_water']
    chl_TW_obs = TW_data['mass_concentration_of_chlorophyll_a_in_sea_water']
    time_TW = TW_data['cast_start_time']

In [5]:

oxy_TW_obs.shape

Out[5]:

(2878, 141)

In [10]:

fig, ax = plt.subplots(5, 1, figsize=(15, 18))
temp_TW_obs.plot(ax=ax[0], y="sea_water_pressure", yincrease=False, cmap='gist_rainbow_r', vmin=3, vmax=21, ylim=(30, 0));
salt_TW_obs.plot(ax=ax[1], y="sea_water_pressure", yincrease=False, cmap= 'gist_rainbow_r', vmin=16, vmax=31, ylim=(30, 0));
oxy_TW_obs.plot(ax=ax[2], y="sea_water_pressure", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=400, ylim=(30, 0));
chl_TW_obs.plot(ax=ax[3], y="sea_water_pressure", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=25, ylim=(30, 0));
no3_TW_obs.plot(ax=ax[4], y="sea_water_pressure", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=50, ylim=(30, 0));

ORCA RAWdata - https://nwem.apl.washington.edu/prod_DataReq.shtml

In [112]:

HC_HP_data = loadmat('/ocean/atall/MOAD/Obs/PugetSound/ORCA_rawdata/HC_HP_CTD_data_bin_web.mat') # Hood Canal - Hoodsport
HC_TW_data = loadmat('/ocean/atall/MOAD/Obs/PugetSound/ORCA_rawdata/HC_TW_CTD_data_bin_web.mat') # Hood Canal - Lynch Cove - Twanoh
HC_DB_data = loadmat('/ocean/atall/MOAD/Obs/PugetSound/ORCA_rawdata/HC_DB_CTD_data_bin_web.mat') # Hood Canal - Duckabush
HC_DK_data = loadmat('/ocean/atall/MOAD/Obs/PugetSound/ORCA_rawdata/HC_DK_CTD_data_bin_web.mat') # Hood Canal - North buoy
HC_NB_data = loadmat('/ocean/atall/MOAD/Obs/PugetSound/ORCA_rawdata/HC_NB_CTD_data_bin_web.mat') # Hood Canal - Hoodsport
CI_data = loadmat('/ocean/atall/MOAD/Obs/PugetSound/ORCA_rawdata/CI_CTD_data_bin_web.mat') # Carr Inlet
PW_data = loadmat('/ocean/atall/MOAD/Obs/PugetSound/ORCA_rawdata/PW_CTD_data_bin_web.mat') # Point Wells


time_HP = HC_HP_data.get('Btime')
time_TW = HC_TW_data.get('Btime')
time_DB = HC_DB_data.get('Btime')
time_DK = HC_DK_data.get('Btime')
time_NB = HC_NB_data.get('Btime')
time_CI = CI_data.get('Btime')
time_PW = PW_data.get('Btime')

depth_HP_obs = HC_HP_data.get('Bdepth')
depth_TW_obs = HC_TW_data.get('Bdepth')
depth_DB_obs = HC_DB_data.get('Bdepth')
depth_DK_obs = HC_DK_data.get('Bdepth')
depth_NB_obs = HC_NB_data.get('Bdepth')
depth_CI_obs = CI_data.get('Bdepth')
depth_PW_obs = PW_data.get('Bdepth')

oxy_HP_obs = HC_HP_data.get('Boxy_umol')
oxy_TW_obs = HC_TW_data.get('Boxy_umol')
oxy_DB_obs = HC_DB_data.get('Boxy_umol')
oxy_DK_obs = HC_DK_data.get('Boxy_umol')
oxy_NB_obs = HC_NB_data.get('Boxy_umol')
oxy_CI_obs = CI_data.get('Boxy_umol')
oxy_PW_obs = PW_data.get('Boxy_umol')

In [70]:

time_HP.shape

Out[70]:

(120, 21358)

In [118]:

depth_TW_obs = np.ma.masked_array(depth_TW_obs)
oxy_TW_obs = np.ma.masked_array(oxy_TW_obs)
#time_TW = np.ma.masked_invalid(time_TW)
#plt.contourf(time_HP,-depth_HP_obs,oxy_HP_obs, levels=[10, 30, 500])
plt.pcolormesh(time_TW)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[118], line 5
      2 oxy_TW_obs = np.ma.masked_array(oxy_TW_obs)
      3 #time_TW = np.ma.masked_invalid(time_TW)
      4 #plt.contourf(time_HP,-depth_HP_obs,oxy_HP_obs, levels=[10, 30, 500])
----> 5 plt.pcolormesh(time_TW, depth_TW_obs, oxy_TW_obs, cmap='RdBu')

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/pyplot.py:2773, in pcolormesh(alpha, norm, cmap, vmin, vmax, shading, antialiased, data, *args, **kwargs)
   2768 @_copy_docstring_and_deprecators(Axes.pcolormesh)
   2769 def pcolormesh(
   2770         *args, alpha=None, norm=None, cmap=None, vmin=None,
   2771         vmax=None, shading=None, antialiased=False, data=None,
   2772         **kwargs):
-> 2773     __ret = gca().pcolormesh(
   2774         *args, alpha=alpha, norm=norm, cmap=cmap, vmin=vmin,
   2775         vmax=vmax, shading=shading, antialiased=antialiased,
   2776         **({"data": data} if data is not None else {}), **kwargs)
   2777     sci(__ret)
   2778     return __ret

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/__init__.py:1442, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1439 @functools.wraps(func)
   1440 def inner(ax, *args, data=None, **kwargs):
   1441     if data is None:
-> 1442         return func(ax, *map(sanitize_sequence, args), **kwargs)
   1444     bound = new_sig.bind(ax, *args, **kwargs)
   1445     auto_label = (bound.arguments.get(label_namer)
   1446                   or bound.kwargs.get(label_namer))

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/axes/_axes.py:6220, in Axes.pcolormesh(self, alpha, norm, cmap, vmin, vmax, shading, antialiased, *args, **kwargs)
   6217 shading = shading.lower()
   6218 kwargs.setdefault('edgecolors', 'none')
-> 6220 X, Y, C, shading = self._pcolorargs('pcolormesh', *args,
   6221                                     shading=shading, kwargs=kwargs)
   6222 coords = np.stack([X, Y], axis=-1)
   6223 # convert to one dimensional array, except for 3D RGB(A) arrays

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/axes/_axes.py:5717, in Axes._pcolorargs(self, funcname, shading, *args, **kwargs)
   5715 if funcname == 'pcolormesh':
   5716     if np.ma.is_masked(X) or np.ma.is_masked(Y):
-> 5717         raise ValueError(
   5718             'x and y arguments to pcolormesh cannot have '
   5719             'non-finite values or be of type '
   5720             'numpy.ma.core.MaskedArray with masked values')
   5721     # safe_masked_invalid() returns an ndarray for dtypes other
   5722     # than floating point.
   5723     if isinstance(X, np.ma.core.MaskedArray):

ValueError: x and y arguments to pcolormesh cannot have non-finite values or be of type numpy.ma.core.MaskedArray with masked values

In [52]:

con_list = [[element for element in upperElement] for upperElement in HC_HP_data['__globals__']]
#con_list.shape

In [55]:

# zip provides us with both the x and y in a tuple.
newData = list(zip(con_list))
columns = ['obj_contour_x', 'obj_contour_y']
df = pd.DataFrame(newData)
print(df)

Empty DataFrame
Columns: []
Index: []

In [51]:

df.plot()

Out[51]:

<Axes: >

In [11]:

with xr.open_dataset('/home/sallen/MEOPAR/grid/mesh_mask202108.nc') as mesh:
    tmask = mesh.tmask
    mbathy = mesh.mbathy

In [12]:

# Constant and data rang
syear = 2022
eyear = 2022
smonth = 8
emonth = 12
startdate = arrow.get(syear, smonth, 1)
enddate = arrow.get(eyear, emonth, 31)
print (startdate, enddate)

2022-08-01T00:00:00+00:00 2022-12-31T00:00:00+00:00

In [13]:

with xr.open_dataset('/results2/SalishSea/nowcast-green.202111/01apr23/SalishSea_1d_20230401_20230401_chem_T.nc') as data_chem:
    oxy = data_chem['dissolved_oxygen']
with xr.open_dataset('/results2/SalishSea/nowcast-green.202111/01apr23/SalishSea_1d_20230401_20230401_biol_T.nc') as data_biol:
    no3 = data_biol['nitrate']
    don = data_biol['dissolved_organic_nitrogen']
    pon = data_biol['particulate_organic_nitrogen']
    dia = data_biol['diatoms']
    fla = data_biol['flagellates']
with xr.open_dataset('/results2/SalishSea/nowcast-green.202111/01apr23/SalishSea_1d_20230401_20230401_grid_T.nc') as data_grid:
    temp = data_grid['votemper']
    salt = data_grid['vosaline']
#    dens = data_grid['sigma_theta'] # almost same as salinity

In [14]:

tmask2_oxy = tmask.assign_coords(z=("z", oxy.deptht.data))
tmask2_oxy = tmask2_oxy.rename({"z": "deptht"})
tmask2_no3 = tmask.assign_coords(z=("z", no3.deptht.data))
tmask2_no3 = tmask2_no3.rename({"z": "deptht"})
tmask2_don = tmask.assign_coords(z=("z", don.deptht.data))
tmask2_don = tmask2_don.rename({"z": "deptht"})
tmask2_pon = tmask.assign_coords(z=("z", pon.deptht.data))
tmask2_pon = tmask2_pon.rename({"z": "deptht"})
tmask2_dia = tmask.assign_coords(z=("z", dia.deptht.data))
tmask2_dia = tmask2_dia.rename({"z": "deptht"})
tmask2_fla = tmask.assign_coords(z=("z", fla.deptht.data))
tmask2_fla = tmask2_fla.rename({"z": "deptht"})
tmask2_temp = tmask.assign_coords(z=("z", temp.deptht.data))
tmask2_temp = tmask2_temp.rename({"z": "deptht"})
tmask2_salt = tmask.assign_coords(z=("z", salt.deptht.data))
tmask2_salt = tmask2_salt.rename({"z": "deptht"})
#tmask2_dens = tmask.assign_coords(z=("z", dens.deptht.data))
#tmask2_dens = tmask2_dens.rename({"z": "deptht"})

In [15]:

print (oxy.deptht)

<xarray.DataArray 'deptht' (deptht: 40)>
array([  0.5     ,   1.500003,   2.500011,   3.500031,   4.500071,   5.500151,
         6.50031 ,   7.500623,   8.501236,   9.502433,  10.504766,  11.509312,
        12.518167,  13.535412,  14.568982,  15.634288,  16.761173,  18.007135,
        19.481785,  21.389978,  24.100256,  28.229916,  34.685757,  44.517723,
        58.484333,  76.58559 ,  98.06296 , 121.866516, 147.08946 , 173.11449 ,
       199.57304 , 226.2603  , 253.06664 , 279.93454 , 306.8342  , 333.75018 ,
       360.67453 , 387.6032  , 414.5341  , 441.4661  ], dtype=float32)
Coordinates:
  * deptht   (deptht) float32 0.5 1.5 2.5 3.5 4.5 ... 360.7 387.6 414.5 441.5
Attributes:
    standard_name:  sea_floor_depth
    long_name:      Sea Floor Depth
    units:          metres
    positive:       down

In [16]:

first = True
i1, i2, j1, j2 = (95, 200, 68, 160)  # Hood Canal all
#i1, i2, j1, j2 = (95, 170, 76, 120)  # Central Hood Canal
#i1, i2, j1, j2 = (106, 119, 69, 79)  # 
#i1lc, i2lc, j1lc, j2lc = (120, 150, 68, 75)  # Lynch Cove
i1hs, i2hs, j1hs, j2hs = (107, 110, 83, 86)           # Hoodsport buoy point in Lynch Cove (Hood Canal) - 47°25'18.5"N 123°06'45.3"W
i1lc, i2lc, j1lc, j2lc = (122, 125, 70, 73)  # Twanoh buoy point in Lynch Cove (Hood Canal) - 47.37°N, 123.01°W
for ii, day in enumerate(arrow.Arrow.range('day', startdate, enddate)):
    year = day.year
    yr2 = day.strftime("%y")
    month = day.month
    Month = day.strftime("%b").lower()
    day = day.day
    # set up filename to follow NEMO conventions
    filename_chem = f'/results2/SalishSea/nowcast-green.202111/{day:02}{Month}{yr2}/SalishSea_1d_{year}{month:02}{day:02}_{year}{month:02}{day:02}_chem_T.nc'
    with xr.open_dataset(filename_chem) as ds_chem:
        oxyhs = ds_chem['dissolved_oxygen'][0, :, j1hs:j2hs, i1hs:i2hs].where(tmask2_oxy[0, :, j1hs:j2hs, i1hs:i2hs] == 1).mean(axis=1).mean(axis=1)
        oxylc = ds_chem['dissolved_oxygen'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_oxy[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)
        if first:
            first = False
            oxyhs_ts = oxyhs
            oxylc_ts = oxylc
        else:
            oxyhs_ts = xr.concat([oxyhs_ts, oxyhs], dim='time_counter')
            oxylc_ts = xr.concat([oxylc_ts, oxylc], dim='time_counter')

In [17]:

oxylc_ts.shape

Out[17]:

(153, 40)

In [18]:

first = True
i1, i2, j1, j2 = (95, 200, 68, 160)  # Hood Canal all
#i1, i2, j1, j2 = (95, 170, 76, 120)  # Central Hood Canal
#i1, i2, j1, j2 = (106, 119, 69, 79)  # 
#i1lc, i2lc, j1lc, j2lc = (120, 150, 68, 75)  # Lynch Cove
i1hs, i2hs, j1hs, j2hs = (107, 110, 83, 86)           # Hoodsport buoy point in Lynch Cove (Hood Canal) - 47°25'18.5"N 123°06'45.3"W
i1lc, i2lc, j1lc, j2lc = (122, 125, 70, 73)  # Twanoh buoy point in Lynch Cove (Hood Canal) - 47.37°N, 123.01°W
for ii, day in enumerate(arrow.Arrow.range('day', startdate, enddate)):
    year = day.year
    yr2 = day.strftime("%y")
    month = day.month
    Month = day.strftime("%b").lower()
    day = day.day
    # set up filename to follow NEMO conventions
    filename_biol = f'/results2/SalishSea/nowcast-green.202111/{day:02}{Month}{yr2}/SalishSea_1d_{year}{month:02}{day:02}_{year}{month:02}{day:02}_biol_T.nc'
    with xr.open_dataset(filename_biol) as ds_biol:
        no3hs = ds_biol['nitrate'][0, :, j1hs:j2hs, i1hs:i2hs].where(tmask2_no3[0, :, j1hs:j2hs, i1hs:i2hs] == 1).mean(axis=1).mean(axis=1)
        no3lc = ds_biol['nitrate'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_no3[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)
        donlc = ds_biol['dissolved_organic_nitrogen'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_don[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)
        donhs = ds_biol['dissolved_organic_nitrogen'][0, :, j1hs:j2hs, i1hs:i2hs].where(tmask2_don[0, :, j1hs:j2hs, i1hs:i2hs] == 1).mean(axis=1).mean(axis=1)
        ponlc = ds_biol['particulate_organic_nitrogen'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_pon[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)
        ponhs = ds_biol['particulate_organic_nitrogen'][0, :, j1hs:j2hs, i1hs:i2hs].where(tmask2_pon[0, :, j1hs:j2hs, i1hs:i2hs] == 1).mean(axis=1).mean(axis=1)
        dialc = ds_biol['diatoms'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_dia[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)
        diahs = ds_biol['diatoms'][0, :, j1hs:j2hs, i1hs:i2hs].where(tmask2_dia[0, :, j1hs:j2hs, i1hs:i2hs] == 1).mean(axis=1).mean(axis=1)
        flalc = ds_biol['flagellates'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_fla[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)
        flahs = ds_biol['flagellates'][0, :, j1hs:j2hs, i1hs:i2hs].where(tmask2_fla[0, :, j1hs:j2hs, i1hs:i2hs] == 1).mean(axis=1).mean(axis=1)
        if first:
            first = False
            no3hs_ts = no3hs
            no3lc_ts = no3lc
            donlc_ts = donlc
            donhs_ts = donhs
            ponlc_ts = ponlc
            ponhs_ts = ponhs
            dialc_ts = dialc
            diahs_ts = diahs
            flalc_ts = flalc
            flahs_ts = flahs
        else:
            no3hs_ts = xr.concat([no3hs_ts, no3hs], dim='time_counter')
            no3lc_ts = xr.concat([no3lc_ts, no3lc], dim='time_counter')
            donlc_ts = xr.concat([donlc_ts, donlc], dim='time_counter')
            donhs_ts = xr.concat([donhs_ts, donhs], dim='time_counter')
            ponlc_ts = xr.concat([ponlc_ts, ponlc], dim='time_counter')
            ponhs_ts = xr.concat([ponhs_ts, ponhs], dim='time_counter')
            dialc_ts = xr.concat([dialc_ts, dialc], dim='time_counter')
            diahs_ts = xr.concat([diahs_ts, diahs], dim='time_counter')
            flalc_ts = xr.concat([flalc_ts, flalc], dim='time_counter')
            diahs_ts = xr.concat([diahs_ts, diahs], dim='time_counter')

In [19]:

chllc_ts = (dialc_ts + flalc_ts) * 2
chlhs_ts = (diahs_ts + flahs_ts) * 2

In [20]:

first = True
i1, i2, j1, j2 = (95, 200, 68, 160)  # Hood Canal all
#i1, i2, j1, j2 = (95, 170, 76, 120)  # Central Hood Canal
#i1, i2, j1, j2 = (106, 119, 69, 79)  # 
#i1lc, i2lc, j1lc, j2lc = (120, 150, 68, 75)  # Lynch Cove
i1hs, i2hs, j1hs, j2hs = (107, 110, 83, 86)           # Hoodsport buoy point in Lynch Cove (Hood Canal) - 47°25'18.5"N 123°06'45.3"W
i1lc, i2lc, j1lc, j2lc = (122, 125, 70, 73)  # Twanoh buoy point in Lynch Cove (Hood Canal) - 47.37°N, 123.01°W
for ii, day in enumerate(arrow.Arrow.range('day', startdate, enddate)):
    year = day.year
    yr2 = day.strftime("%y")
    month = day.month
    Month = day.strftime("%b").lower()
    day = day.day
    # set up filename to follow NEMO conventions
    filename_grid = f'/results2/SalishSea/nowcast-green.202111/{day:02}{Month}{yr2}/SalishSea_1d_{year}{month:02}{day:02}_{year}{month:02}{day:02}_grid_T.nc'
    with xr.open_dataset(filename_grid) as ds_grid:
        temphs = ds_grid['votemper'][0, :, j1hs:j2hs, i1hs:i2hs].where(tmask2_temp[0, :, j1hs:j2hs, i1hs:i2hs] == 1).mean(axis=1).mean(axis=1)
        templc = ds_grid['votemper'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_temp[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)
        salths = ds_grid['vosaline'][0, :, j1hs:j2hs, i1hs:i2hs].where(tmask2_salt[0, :, j1hs:j2hs, i1hs:i2hs] == 1).mean(axis=1).mean(axis=1)
        saltlc = ds_grid['vosaline'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_salt[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)
        #dens = ds_grid['sigma_theta'][0, :, j1:j2, i1:i2].where(tmask2_dens[0, :, j1:j2, i1:i2] == 1).mean(axis=1).mean(axis=1)
        #denslc = ds_grid['sigma_theta'][0, :, j1lc:j2lc, i1lc:i2lc].where(tmask2_dens[0, :, j1lc:j2lc, i1lc:i2lc] == 1).mean(axis=1).mean(axis=1)        
        if first:
            first = False
            temphs_ts = temphs
            templc_ts = templc
            salths_ts = salths
            saltlc_ts = saltlc
            #dens_ts = dens
            #denslc_ts = denslc
        else:
            temphs_ts = xr.concat([temphs_ts, temphs], dim='time_counter')
            templc_ts = xr.concat([templc_ts, templc], dim='time_counter')
            salths_ts = xr.concat([salths_ts, salths], dim='time_counter')
            saltlc_ts = xr.concat([saltlc_ts, saltlc], dim='time_counter')
            #dens_ts = xr.concat([dens_ts, dens], dim='time_counter')
            #denslc_ts = xr.concat([denslc_ts, denslc], dim='time_counter')

In [21]:

fig, ax = plt.subplots(5, 1, figsize=(15, 18))
templc_ts[:, :24].plot(ax=ax[0], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=3, vmax=21, ylim=(30, 0));
saltlc_ts[:, :24].plot(ax=ax[1], y="deptht", yincrease=False, cmap= 'gist_rainbow_r', vmin=16, vmax=31, ylim=(30, 0));
oxylc_ts[:, :24].plot(ax=ax[2], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=400, ylim=(30, 0));
chllc_ts[:, :24].plot(ax=ax[3], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=25, ylim=(30, 0));
no3lc_ts[:, :24].plot(ax=ax[4], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=50, ylim=(30, 0));

ax[0].set_ylabel('Depth (m)')
ax[1].set_ylabel('Depth (m)')
ax[2].set_ylabel('Depth (m)')
ax[3].set_ylabel('Depth (m)')
ax[0].set_xlabel('')
ax[1].set_xlabel('')
ax[2].set_xlabel('')
ax[3].set_xlabel('')
ax[4].set_xlabel('Hood Canal (Twanoh - Lynch Cove)')

ax[0].set_title('Temperature (°C)')
ax[1].set_title('Salinity (g/kg)')
ax[2].set_title('Oxygen (µmol/L)')
ax[3].set_title('Chl (mg/m3)')
ax[4].set_title('Nitrate (µmol/L)')

# Make ticks on occurrences of each month:
ax[0].xaxis.set_major_locator(mdates.MonthLocator())
ax[1].xaxis.set_major_locator(mdates.MonthLocator())
ax[2].xaxis.set_major_locator(mdates.MonthLocator())
ax[3].xaxis.set_major_locator(mdates.MonthLocator())
ax[4].xaxis.set_major_locator(mdates.MonthLocator())
# Get only the month to show in the x-axis:
ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[1].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[2].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[3].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[4].xaxis.set_major_formatter(mdates.DateFormatter('%b'))

In [22]:

fig, ax = plt.subplots(5, 1, figsize=(15, 18))
temphs_ts[:, :24].plot(ax=ax[0], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=3, vmax=21, ylim=(30, 0));
salths_ts[:, :24].plot(ax=ax[1], y="deptht", yincrease=False, cmap= 'gist_rainbow_r', vmin=16, vmax=31, ylim=(30, 0));
oxyhs_ts[:, :24].plot(ax=ax[2], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=400, ylim=(30, 0));
chlhs_ts[:, :24].plot(ax=ax[3], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=25, ylim=(30, 0));
no3hs_ts[:, :24].plot(ax=ax[4], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=50, ylim=(30, 0));

ax[0].set_ylabel('Depth (m)')
ax[1].set_ylabel('Depth (m)')
ax[2].set_ylabel('Depth (m)')
ax[3].set_ylabel('Depth (m)')
ax[0].set_xlabel('')
ax[1].set_xlabel('')
ax[2].set_xlabel('')
ax[3].set_xlabel('')
ax[4].set_xlabel('Hood Canal (Hoodsport - Central Hood Canal)')

ax[0].set_title('Temperature (°C)')
ax[1].set_title('Salinity (g/kg)')
ax[2].set_title('Oxygen (µmol/L)')
ax[3].set_title('Chl (mg/m3)')
ax[4].set_title('Nitrate (µmol/L)')

# Make ticks on occurrences of each month:
ax[0].xaxis.set_major_locator(mdates.MonthLocator())
ax[1].xaxis.set_major_locator(mdates.MonthLocator())
ax[2].xaxis.set_major_locator(mdates.MonthLocator())
ax[3].xaxis.set_major_locator(mdates.MonthLocator())
ax[4].xaxis.set_major_locator(mdates.MonthLocator())
# Get only the month to show in the x-axis:
ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[1].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[2].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[3].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[4].xaxis.set_major_formatter(mdates.DateFormatter('%b'))

In [23]:

fig, ax = plt.subplots(3, 1, figsize=(18, 10))
donlc_ts[:, :22].plot(ax=ax[0], y="deptht", yincrease=False, cmap= 'gist_rainbow_r', vmin=0, vmax=5, ylim=(30, 0));
ponlc_ts[:, :22].plot(ax=ax[1], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=2, ylim=(30, 0));
oxylc_ts[:, :22].plot(ax=ax[2], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=60, vmax=300, ylim=(30, 0));
ax[2].set_xlim(startdate, enddate)

ax[0].set_ylabel('Depth (m)')
ax[1].set_ylabel('Depth (m)')
ax[2].set_ylabel('Depth (m)')
ax[0].set_xlabel('')
ax[1].set_xlabel('')
ax[2].set_xlabel('Hood Canal (Twanoh - Lynch Cove)')

ax[0].set_title('DON (µmol/L)')
ax[1].set_title('PON (µmol/L)')
ax[2].set_title('Oxygen (µmol/L)')

# Make ticks on occurrences of each month:
ax[0].xaxis.set_major_locator(mdates.MonthLocator())
ax[1].xaxis.set_major_locator(mdates.MonthLocator())
ax[2].xaxis.set_major_locator(mdates.MonthLocator())
# Get only the month to show in the x-axis:
ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[1].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[2].xaxis.set_major_formatter(mdates.DateFormatter('%b'))

In [24]:

fig, ax = plt.subplots(3, 1, figsize=(18, 10))
donhs_ts[:, :22].plot(ax=ax[0], y="deptht", yincrease=False, cmap= 'gist_rainbow_r', vmin=0, vmax=5, ylim=(30, 0));
ponhs_ts[:, :22].plot(ax=ax[1], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=0, vmax=2, ylim=(30, 0));
oxyhs_ts[:, :22].plot(ax=ax[2], y="deptht", yincrease=False, cmap='gist_rainbow_r', vmin=60, vmax=300, ylim=(30, 0));
ax[2].set_xlim(startdate, enddate)

ax[0].set_ylabel('Depth (m)')
ax[1].set_ylabel('Depth (m)')
ax[2].set_ylabel('Depth (m)')
ax[0].set_xlabel('')
ax[1].set_xlabel('')
ax[2].set_xlabel('Hood Canal (Twanoh - Lynch Cove)')

ax[0].set_title('DON (µmol/L)')
ax[1].set_title('PON (µmol/L)')
ax[2].set_title('Oxygen (µmol/L)')

# Make ticks on occurrences of each month:
ax[0].xaxis.set_major_locator(mdates.MonthLocator())
ax[1].xaxis.set_major_locator(mdates.MonthLocator())
ax[2].xaxis.set_major_locator(mdates.MonthLocator())
# Get only the month to show in the x-axis:
ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[1].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax[2].xaxis.set_major_formatter(mdates.DateFormatter('%b'))

In [25]:

fig, ax = plt.subplots(1, 2, figsize=(10, 9))
ax[0].plot(no3hs_ts[91:273, :22], oxyhs_ts[91:273, :22], '.');
ax[1].plot(no3lc_ts[91:273, :22], oxylc_ts[91:273, :22], '.');
ax[0].set_xlim(-1, 30)
ax[1].set_xlim(-1, 30)
ax[0].set_ylim(90, 380)
ax[1].set_ylim(90, 380)
ax[0].set_ylabel('Dissolved_oxygen')
ax[0].set_xlabel('NO3')
ax[1].set_xlabel('NO3')
ax[0].text(10, 350, 'Hoodsport (0-30 m, Apr-Sep)')
ax[1].text(10, 350, 'Lynch Cove (0-30 m, Apr-Sep)')

Out[25]:

Text(10, 350, 'Lynch Cove (0-30 m, Apr-Sep)')

In [26]:

fig, ax = plt.subplots(6, 1, figsize=(15, 14))
#specify x-axis locations
#x_ticks = [1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335]
x_ticks = [1, 366, 732, 1097, 1462, 1827, 2193, 2558]
#specify x-axis labels
#x_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] 
x_labels = ['2007', '2008', '2009', '2010', '2011', '2012', '2013', '2014'] 
t1 = np.mean(templc_ts[:, :4], axis=1)-np.mean(templc_ts[:, :4], axis=(0, 1))
t2 = np.mean(templc_ts[:, 5:], axis=1)-np.mean(templc_ts[:, 5:], axis=(0, 1))
plt.setp(ax, xticks=x_ticks, xticklabels=x_labels)
ax[0].plot(t1,'b')
#ax[0].set_title('Top 5m')
ax[0].set_ylabel('Temp (°C)')
ax[0].plot(t2,'g')
ax[0].grid()
s1 = np.mean(saltlc_ts[:, :4], axis=1)-np.mean(saltlc_ts[:, :4], axis=(0, 1))
s2 = np.mean(saltlc_ts[:, 5:], axis=1)-np.mean(saltlc_ts[:, 5:], axis=(0, 1))
ax[1].plot(s1,'b')
ax[1].set_ylabel('Salt (g/kg)')
ax[1].plot(s2,'g')
ax[1].grid()
n1 = np.mean(no3lc_ts[:, :4], axis=1)-np.mean(no3lc_ts[:, :4], axis=(0, 1))
n2 = np.mean(no3lc_ts[:, 5:], axis=1)-np.mean(no3lc_ts[:, 5:], axis=(0, 1))
ax[2].plot(n1,'b')
ax[2].set_ylabel('NO3 (µmol/L)')
ax[2].plot(n2,'g')
ax[2].grid()
d1 = np.mean(donlc_ts[:, :4], axis=1)-np.mean(donlc_ts[:, :4], axis=(0, 1))
d2 = np.mean(donlc_ts[:, 5:], axis=1)-np.mean(donlc_ts[:, 5:], axis=(0, 1))
ax[3].plot(d1,'b')
ax[3].set_ylabel('DON (µmol/L)')
ax[3].plot(d2,'g')
ax[3].grid()
p1 = np.mean(ponlc_ts[:, :4], axis=1)-np.mean(ponlc_ts[:, :4], axis=(0, 1))
p2 = np.mean(ponlc_ts[:, 5:], axis=1)-np.mean(ponlc_ts[:, 5:], axis=(0, 1))
ax[4].plot(p1,'b')
ax[4].set_ylabel('PON (µmol/L)')
ax[4].plot(p2,'g')
ax[4].grid()
o1 = np.mean(oxylc_ts[:, :4], axis=1)-np.mean(oxylc_ts[:, :4], axis=(0, 1))
o2 = np.mean(oxylc_ts[:, 5:], axis=1)-np.mean(oxylc_ts[:, 5:], axis=(0, 1))
ax[5].plot(o1,'b')
ax[5].set_ylabel('Oxy (µmol/L)')
ax[5].plot(o2,'g')
ax[5].grid()

In [27]:

fig, ax = plt.subplots(6, 1, figsize=(15, 14))
#specify x-axis locations
#x_ticks = [1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335]
x_ticks = [1, 366, 732, 1097, 1462, 1827, 2193, 2558]
#specify x-axis labels
#x_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] 
x_labels = ['2007', '2008', '2009', '2010', '2011', '2012', '2013', '2014'] 
t1 = np.mean(templc_ts[:, :4], axis=1)-np.mean(templc_ts[:, :4], axis=(0, 1))
t2 = np.mean(templc_ts[:, 5:], axis=1)-np.mean(templc_ts[:, 5:], axis=(0, 1))
plt.setp(ax, xticks=x_ticks, xticklabels=x_labels)
ax[0].set_title('Below 5m')
ax[0].set_ylabel('Temp (°C)')
ax[0].plot(t2,'g')
ax[0].grid()
s1 = np.mean(saltlc_ts[:, :4], axis=1)-np.mean(saltlc_ts[:, :4], axis=(0, 1))
s2 = np.mean(saltlc_ts[:, 5:], axis=1)-np.mean(saltlc_ts[:, 5:], axis=(0, 1))
ax[1].set_ylabel('Salt (g/kg)')
ax[1].plot(s2,'g')
ax[1].grid()
n1 = np.mean(no3lc_ts[:, :4], axis=1)-np.mean(no3lc_ts[:, :4], axis=(0, 1))
n2 = np.mean(no3lc_ts[:, 5:], axis=1)-np.mean(no3lc_ts[:, 5:], axis=(0, 1))
ax[2].set_ylabel('NO3 (µmol/L)')
ax[2].plot(n2,'g')
ax[2].grid()
d1 = np.mean(donlc_ts[:, :4], axis=1)-np.mean(donlc_ts[:, :4], axis=(0, 1))
d2 = np.mean(donlc_ts[:, 5:], axis=1)-np.mean(donlc_ts[:, 5:], axis=(0, 1))
ax[3].set_ylabel('DON (µmol/L)')
ax[3].plot(d2,'g')
ax[3].grid()
p1 = np.mean(ponlc_ts[:, :4], axis=1)-np.mean(ponlc_ts[:, :4], axis=(0, 1))
p2 = np.mean(ponlc_ts[:, 5:], axis=1)-np.mean(ponlc_ts[:, 5:], axis=(0, 1))
ax[4].set_ylabel('PON (µmol/L)')
ax[4].plot(p2,'g')
ax[4].grid()
o1 = np.mean(oxylc_ts[:, :4], axis=1)-np.mean(oxylc_ts[:, :4], axis=(0, 1))
o2 = np.mean(oxylc_ts[:, 5:], axis=1)-np.mean(oxylc_ts[:, 5:], axis=(0, 1))
ax[5].set_ylabel('Oxy (µmol/L)')
ax[5].plot(o2,'g')
ax[5].grid()

In [28]:

i1, i2, j1, j2 = (95, 200, 68, 160)  # Hood Canal all
fig, ax = plt.subplots(1, 1, figsize=(5, 9))
oxy[0, mbathy[0, :, :]-1, :, :].where(tmask[0, mbathy[0, :, :]-1, :, :] == 1).plot(ax=ax, cmap=cm.oxy, vmin=20, vmax=250)
viz_tools.set_aspect(ax);
# Oxygen in the Hood Canal (j=60:120, i=100:160)
left, bottom, width, height = (95, 68, 105, 100)
rect1=mpatches.Rectangle((left,bottom),width,height, 
                        fill=False,
                        #alpha=0.1
                        color="purple",
                       linewidth=2)
                       #facecolor="red")
left, bottom, width, height = (180, 330, 35, 60)
rect2=mpatches.Rectangle((left,bottom),width,height, 
                        fill=False,
                        #alpha=0.1
                        color="purple",
                       linewidth=2)
                       #facecolor="red")
plt.gca().add_patch(rect1)
plt.gca().add_patch(rect2)

fig.suptitle('Bottom Oxygen - 202111-nowcast - Year 2011');

#120, 150, 68, 75  # Lynch Cove mooring
#107, 110, 83, 86  # Hoodsport mooring
ax.set_ylim(75, 100)
ax.set_xlim(100, 120)

Out[28]:

(100.0, 120.0)

In [ ]: