<jupyter_start><jupyter_text>202007Jp2 bloom timing<jupyter_code>import datetime as dt
from erddapy import ERDDAP
import glob
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import pandas as pd
from pathlib import Path
import xarray as xr

from salishsea_tools import places

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Set up Model Results Getting<jupyter_code>directory = '202007J-p2'
dir_path = Path(f"/data/sallen/results/MEOPAR/202007/{directory}")
drop_vars = [
    "time_centered",
    "time_centered_bounds",
    "time_counter_bounds",
    "nav_lat",
    "bounds_lat",
    "nav_lon",
    "bounds_lon",
    "area",
    "deptht_bounds",
]
chunks = {
    "time_counter": 1,
    "deptht": 1,
}
def prep_timeseries(dir_path, file_pattern, n_days): 
    ds_files = sorted(dir_path.glob(file_pattern))
    ds = xr.open_mfdataset(
        ds_files[:n_days],
        chunks=chunks,
        drop_variables=drop_vars,
    )
    S3y, S3x = places.PLACES['S3']['NEMO grid ji'][0], places.PLACES['S3']['NEMO grid ji'][1]
    fourkm = 8
    timeseries = ds.isel(deptht = 1, y = slice(S3y-fourkm, S3y+fourkm+1), 
                             x =  slice(S3x-fourkm, S3x+fourkm+1))
    
    return timeseries<jupyter_output><empty_output><jupyter_text>Set up Observations Getting<jupyter_code>server = "https://salishsea.eos.ubc.ca/erddap"

protocol = "tabledap"

dataset_id = "ubcONCTWDP1mV18-01"

response = "nc"

variables = [
    "latitude",
    "longitude",
    "chlorophyll",
    "o2_concentration_corrected",
    "time",
]

variables2 = [
    "latitude",
    "longitude",
    "o2_concentration_corrected",
    "time",
]

    

fourkmlat = 4/110.574
fourkmlon = 4/(111.320*np.cos(50*np.pi/180.))
lon, lat = places.PLACES['S3']['lon lat']<jupyter_output><empty_output><jupyter_text>Plotting Functions<jupyter_code>def model_plot(timeseries2015):
    fig, axs = plt.subplots(1, 2, figsize=(15, 4))
    for i in range(timeseries2015.x.shape[0]):
        for j in range(timeseries2015.y.shape[0]):
            timeseries2015.diatoms[:, j, i].plot(ax=axs[0], color='blue', marker='o')
            timeseries2015.diatoms[:, j, i].plot(ax=axs[1], color='blue', marker='o')
    ax2 = axs[0].twinx()
    for i in range(timeseries2015.x.shape[0]):
        for j in range(timeseries2015.y.shape[0]):
            timeseries2015.nitrate[:, j, i].plot(ax=ax2, color='lightblue', marker='o')
    axs[0].grid()
    for i in range(timeseries2015.x.shape[0]):
        for j in range(timeseries2015.y.shape[0]):
            (timeseries2015.diatoms+timeseries2015.flagellates+timeseries2015.ciliates)[:, j, i].plot(ax=axs[1], color='purple', marker='o', alpha=0.5);
def model_obs_comparison(timeseries2015, obs_pd, N2Chl=2.5, oxy=False):
    fig, axs = plt.subplots(1, 1, figsize=(7, 4))
    if oxy:
        ax2 = axs
        ax = ax2.twinx()
        obs_pd['o2_concentration_corrected (ml/l)'].plot(ax=ax2, marker='o', color='lime', linewidth=0, alpha=0.4, zorder=6);
    else:
        ax = axs
        obs_pd['chlorophyll (ug/l)'].plot(ax=ax, marker='o', color='green', linewidth=0, zorder=2)
    for i in range(timeseries2015.x.shape[0]):
        for j in range(timeseries2015.y.shape[0]):
            ((timeseries2015.diatoms+timeseries2015.flagellates+timeseries2015.ciliates)[:, j, i]*N2Chl).plot(ax=ax, color='purple', 
                                                                                                              marker='.', alpha=0.5, linewidth=0, zorder=3);<jupyter_output><empty_output><jupyter_text>2015  Model Results<jupyter_code>year = 2015
months_regex = '0[23]'
file_pattern = f"SalishSea_1d_{year}*ptrc_T*{year}{months_regex}??.nc"
timeseries = prep_timeseries(dir_path, file_pattern, -1)
timeseries.load(scheduler='processes', num_workers=4);
model_plot(timeseries)
timeseries.diatoms[28-1:28-1+10].mean(axis=1).mean(axis=1).plot();
plt.grid()<jupyter_output><empty_output><jupyter_text>Bloom is March 7.  This bloom is two days early.  Observations<jupyter_code>constraints = {
    "time>=": "2015-02-01T00:00:00Z",
    "time<=": "2015-04-01T00:00:00Z",
    "latitude>=": lat - fourkmlat,
    "latitude<=": lat + fourkmlat,
    "longitude>=": lon - fourkmlon,
    "longitude<=": lon + fourkmlon,
}
obs = ERDDAP(server=server, protocol=protocol,)

obs.dataset_id = dataset_id
obs.variables = variables
obs.constraints = constraints
obs_pd = obs.to_pandas(index_col="time (UTC)", parse_dates=True,)
obs2 = ERDDAP(server=server, protocol=protocol,)

obs2.dataset_id = dataset_id
obs2.variables = variables2
obs2.constraints = constraints
obs2_pd = obs2.to_pandas(index_col="time (UTC)", parse_dates=True,).dropna()
model_obs_comparison(timeseries, obs_pd, N2Chl=3)
model_obs_comparison(timeseries, obs_pd, N2Chl=3, oxy=True)
fig, ax = plt.subplots(1, 1, figsize=(10,4))
ax2 = ax.twinx()
obs_pd['o2_concentration_corrected (ml/l)'].groupby(pd.Grouper(freq='1D')).mean().plot(ax=ax)
obs_pd['chlorophyll (ug/l)'].groupby(pd.Grouper(freq='1D')).mean().plot(ax=ax2, color='orange');
ax.grid()<jupyter_output><empty_output><jupyter_text>Bloom is March 9.  2016  Model Results<jupyter_code>year = 2016
months_regex = '0[34]'
file_pattern = f"SalishSea_1d_{year}*ptrc_T*{year}{months_regex}??.nc"
timeseries = prep_timeseries(dir_path, file_pattern, -1)
timeseries.load(scheduler='processes', num_workers=4);
model_plot(timeseries)
timeseries.diatoms[29-1:32-1].mean(axis=1).mean(axis=1)
timeseries.diatoms[28-5:28-5+10].mean(axis=1).mean(axis=1).plot();
plt.grid()<jupyter_output><empty_output><jupyter_text>Bloom date is March 30.  This value is only 9 days late.  Observations<jupyter_code>constraints = {
    "time>=": "2016-03-01T00:00:00Z",
    "time<=": "2016-05-01T00:00:00Z",
    "latitude>=": lat - fourkmlat,
    "latitude<=": lat + fourkmlat,
    "longitude>=": lon - fourkmlon,
    "longitude<=": lon + fourkmlon,
}
obs = ERDDAP(server=server, protocol=protocol,)

obs.dataset_id = dataset_id
obs.variables = variables
obs.constraints = constraints
obs_pd = obs.to_pandas(index_col="time (UTC)", parse_dates=True,).dropna()
obs2 = ERDDAP(server=server, protocol=protocol,)

obs2.dataset_id = dataset_id
obs2.variables = variables2
obs2.constraints = constraints
obs2_pd = obs2.to_pandas(index_col="time (UTC)", parse_dates=True,).dropna()
model_obs_comparison(timeseries, obs_pd, N2Chl=1.4)
model_obs_comparison(timeseries, obs_pd, N2Chl=1.4, oxy=True)
fig, ax = plt.subplots(1, 1, figsize=(10,4))
ax2 = ax.twinx()
obs_pd['chlorophyll (ug/l)'].plot(ax=ax)
obs_pd['o2_concentration_corrected (ml/l)'].plot(ax=ax2, color='tab:orange')
fig, ax = plt.subplots(1, 1, figsize=(10,4))
ax2 = ax.twinx()
obs_pd['chlorophyll (ug/l)'].groupby(pd.Grouper(freq='1D')).mean().plot(ax=ax)
obs_pd['o2_concentration_corrected (ml/l)'].groupby(pd.Grouper(freq='1D')).mean().plot(ax=ax2, color='tab:orange')
obs_pd['chlorophyll (ug/l)'].groupby(pd.Grouper(freq='1D')).mean()[10:29]<jupyter_output><empty_output><jupyter_text>Bloom date is March 21.  2017  Model Results<jupyter_code>year = 2017
months_regex = '0[234]'
file_pattern = f"SalishSea_1d_{year}*ptrc_T*{year}{months_regex}??.nc"
timeseries = prep_timeseries(dir_path, file_pattern, -1)
timeseries.load(scheduler='processes', num_workers=4);
model_plot(timeseries)
timeseries.diatoms[32-1:42-1].mean(axis=1).mean(axis=1)<jupyter_output><empty_output><jupyter_text>Split bloom.  Observations<jupyter_code>constraints = {
    "time>=": "2017-02-01T00:00:00Z",
    "time<=": "2017-05-01T00:00:00Z",
    "latitude>=": lat - fourkmlat,
    "latitude<=": lat + fourkmlat,
    "longitude>=": lon - fourkmlon,
    "longitude<=": lon + fourkmlon,
}
obs = ERDDAP(server=server, protocol=protocol,)

obs.dataset_id = dataset_id
obs.variables = variables
obs.constraints = constraints
obs_pd = obs.to_pandas(index_col="time (UTC)", parse_dates=True,).dropna()
model_obs_comparison(timeseries, obs_pd, N2Chl=3)
obs2 = ERDDAP(server=server, protocol=protocol,)

obs2.dataset_id = dataset_id
obs2.variables = variables2
obs2.constraints = constraints
obs2_pd = obs2.to_pandas(index_col="time (UTC)", parse_dates=True,)
model_obs_comparison(timeseries, obs_pd, N2Chl=3, oxy=True)
fig, ax = plt.subplots(1, 1, figsize=(10,4))
ax2 = ax.twinx()
obs_pd['chlorophyll (ug/l)'].groupby(pd.Grouper(freq='1D')).mean().plot(ax=ax)
obs_pd['o2_concentration_corrected (ml/l)'].groupby(pd.Grouper(freq='1D')).mean().plot(ax=ax2, color='tab:orange');
obs_pd['chlorophyll (ug/l)'].groupby(pd.Grouper(freq='1D')).mean()[30:50]
#Bloom is March 21
timeseries.diatoms[31:90].mean(axis=1).mean(axis=1).plot();
plt.grid()
#Bloom is April 4, this is 14 days late.<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16655108645399908, 'fraction_numerical': 0.0518955154877485, 'mean_word_length': 5.149964463397299, 'pattern_counts': {'":': 20, '<': 50, '<?xml version=': 0, '>': 50, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '313182', 'license': 'Apache-2.0', 'n_tokens_mistral': 3692, 'n_tokens_neox': 3414, 'n_words': 526, 'path': 'notebooks/Tuning/Bloom_Timing_202007Jp2.ipynb', 'repo_name': 'SalishSeaCast/analysis-susan', 'size': '8652'}