<jupyter_start><jupyter_text>Process Sentinel3 Chlorophyll<jupyter_code>imagery_root = '/tmp/nepacs3/'
image_files = glob(os.path.join(imagery_root, "*/*nn*.tif"))
def get_date_from_file(file):
    return(pd.to_datetime(os.path.splitext(file)[0].split('-')[1]))
dates = list(map(get_date_from_file, image_files))
xr_ds = [
    rio.open(image_file) for image_file in image_files
]
len(xr_ds)<jupyter_output><empty_output><jupyter_text>**Crop** each one to a buffered zone 'round the chosen isolate (id:240, the nepac one)<jupyter_code>plankton = pd.read_csv(
    '../data/Phytoplankton_temperature_growth_rate_dataset_2016_01_29/traits_derived_2016_01_29.csv',
    engine='python'
)
isolate_code = 240
isolate = plankton[plankton['isolate.code'] == isolate_code]
isolate
project_merc = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'), # source coordinate system
    pyproj.Proj(init='epsg:3857')) # mercator, to compute area buffer

project_wgs = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:3857'), # source coordinate system
    pyproj.Proj(init='epsg:4326')) # mercator, to compute area buffer


isolate_point_wgs = geometry.Point(
    isolate['isolation.longitude'],
    isolate['isolation.latitude']
)
isolate_point_meters = transform(project, isolate_point_wgs)
merc_500m_buffer = isolate_point_meters.buffer(20000)
clipped_data = []
for ds in xr_ds:
    try:
        clipped_data.append(rio.mask.mask(ds, [merc_500m_buffer], crop=True))
    except Exception as e:
        print(e)
        clipped_data.append(e)
    
    
len(clipped_data)<jupyter_output><empty_output><jupyter_text>Turn them all to masked arrays<jupyter_code>ma_clipped = []
for data in clipped_data:
    print(type(data))
    if isinstance(data, Exception): 
        ma_clipped.append(data)
        continue
        
    _ma = np.ma.MaskedArray(
        data = data[0][0], 
        mask = (data[0][1] == 0)
    )
    ma_clipped.append(_ma)<jupyter_output><class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'>
<class 'tuple'><jupyter_text>Compute Mean chl-$\alpha$<jupyter_code>chl = []
for ma in ma_clipped:
    print(isinstance(ma, Exception))
    if isinstance(ma, Exception): 
        chl.append(np.nan)
        continue
    chl.append(np.mean(ma))
    
chl_df = pd.DataFrame({
    'date': dates,
    'chl': chl
})
chl_df.loc[chl_df.chl.apply(lambda x: isinstance(x, np.ma.core.MaskedConstant)), 'chl'] = np.nan
chl_df=chl_df.sort_index()
chl_df.index
plt.scatter(chl_df.date, chl_df.chl)
plt.xticks(rotation=30, ha='right')

len(chl_df[~chl_df.chl.isna()])<jupyter_output><empty_output><jupyter_text>Can we look at performance around this time?<jupyter_code>isolate_performance_perf = xr.open_zarr("../data/isolate_mhw_performance.zarr/isolate_mhw_performance.zarr/")
mhw_performance = isolate_performance_perf.sel(
    isolate=isolate_code, 
    time=slice(chl_df.date.min(), chl_df.date.max())
)
mhw_performance.time.min()
mhw_performance.performance.plot()
plt.ylim([0, 1])
ax2 = plt.twinx()
ax2.scatter(chl_df.date, chl_df.chl, color='red')
plt.xticks(rotation=30, ha='right')
chl_df['date_delta'] = (chl_df.reset_index()['date'] - chl_df.reset_index()['date'].min())  / np.timedelta64(1,'D')
fg = sns.lmplot(x='date_delta', y='chl', data=chl_df.reset_index().dropna(), )
sns.scatterplot(mhw_performance.perf)
fg.axes[0][0]
chl_df[''(chl_df.reset_index()['date'] - chl_df.reset_index()['date'].min())  / np.timedelta64(1,'D')
def lowess_smooth(data):
    from statsmodels.nonparametric.smoothers_lowess import lowess
    
    endog = data.values
    exog = data.index.values
    
    smooth = lowess(endog, exog)
    index, xformed = np.transpose(smooth)

    return pd.Series(xformed, index=pd.to_datetime(index)) 
plt.plot(lowess_smooth(chl_df.dropna(subset=['chl']).set_index('date').chl))
plt.figure(dpi=150)
sns.set(style='ticks')
ax1 = plt.axes()
_p0 = mhw_performance.performance.plot(label='Isolate Performance', ax=ax1)
ax1.set_ylim([0, 1])
ax2 = ax1.twinx()
_p1 = ax2.scatter(chl_df.date, chl_df.chl, color='none', alpha=1, edgecolor='red', label=r'Sentinel-3 Chl-$\alpha$')
plt.xticks(rotation=30, ha='right')
ax2.set_ylabel(r"Chl-$\alpha$ (oc4me)")
_p2 = plt.plot(
    lowess_smooth(chl_df.dropna(subset=['chl']).set_index('date').chl),
    linewidth=2, color='red', linestyle='--',
    label=r'Sentinel-3 Chl-$\alpha$ LOWESS'
)
    
lines, labels = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2, loc='lower right')
sns.despine(right=False)
plt.savefig("nepac_s3_validate.pdf", bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Create Paired plot<jupyter_code>fs = gcsfs.GCSFileSystem(project=GCP_PROJECT_ID, token="../gc-<PASSWORD>")
oisst = xr.open_zarr(fs.get_mapper(OISST_GCP))
oisst = oisst.assign_coords(lon=(((oisst.lon + 180) % 360) - 180)).sortby('lon')
selected_sst = oisst.sel(
    lat = mhw_performance.lat.values.item(),
    lon = mhw_performance.lon.values.item(), 
    time = mhw_performance.time
)
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5), dpi=150)

sst_ax = axes[0]

selected_sst.sst.plot(ax=sst_ax, color='black')
mhw_performance.clim_seas.plot(ax=sst_ax, color='blue')
mhw_performance.clim_thresh.plot(ax=sst_ax, color='green')
sst_ax.set_ylabel("SST [deg C]")
sst_ax.set_title("")
sst_ax.set_title("A) Marine Heat Wave", loc='left')

chl_ax = axes[1]
_p0 = mhw_performance.performance.plot(label='Isolate Performance', ax=chl_ax)
chl_ax.set_ylim([0, 1])
chl_twinax = chl_ax.twinx()
_p1 = chl_twinax.scatter(chl_df.date, chl_df.chl, color='none', alpha=1, edgecolor='red', zorder=-1,  label=r'Sentinel-3 Chl-$\alpha$')
chl_ax.set_xticklabels(chl_ax.get_xticks(), rotation=30, ha='right')
chl_twinax.set_ylabel(r"Chl-$\alpha$ (oc4me)")
_p2 = chl_twinax.plot(
    lowess_smooth(chl_df.dropna(subset=['chl']).set_index('date').chl),
    linewidth=2, color='red', linestyle='--',
    label=r'Sentinel-3 Chl-$\alpha$ LOWESS'
)
    
lines, labels = chl_ax.get_legend_handles_labels()
lines2, labels2 = chl_twinax.get_legend_handles_labels()
chl_twinax.legend(lines + lines2, labels + labels2, loc='center right', fancybox=True, bbox_to_anchor=(1.75, 0.5))

chl_ax.set_title("")
chl_ax.set_title("B) MHW Performance + Validation", loc='left')

sns.despine(right=False)
plt.savefig("nepac_s3_combined.pdf", bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17129564359227406, 'fraction_numerical': 0.02072465811363316, 'mean_word_length': 5.788516746411483, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 73, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '511750', 'license': 'MIT', 'n_tokens_mistral': 3011, 'n_tokens_neox': 2857, 'n_words': 450, 'path': 'validation/nwpacific-example-sentinel.ipynb', 'repo_name': 'HuckleyLab/phyto-mhw', 'size': '7094'}