Commit 0e5748b6 authored by Clément Haëck's avatar Clément Haëck
Browse files

Up method to find bloom peak

add filtering
parent 2e081c51
#!/usr/bin/env ipython
import numpy as np
from os import path
import matplotlib.pyplot as plt
import tol_colors as tc
import xarray as xr
from scipy import signal
import lib
import lib.data.hists
import lib.data.hists as lh
import Plots.util as plot_util
def plot(m, tmax, i=-1):
fig, ax = plt.subplots(figsize=(10, 5))
fig.subplots_adjust(left=0.05, bottom=0.05, right=0.98, top=0.98)
PLOT = True
DAY_RES = 1
DATE_START = '0101'
DATE_STOP = '0830'
# DATE_START = '0901'
# DATE_STOP = '1231'
ZONE = 'GS3_N'
METHOD = 'filter', 15
args = dict(mask=True, kind='2thr')
fixes = dict(
number=2,
scale=30,
coef=0,
thr_lo=5.,
thr_hi=10.,
# Y=list(range(2003, 2015, 3))
)
ds = lh.get_data(args, fixes=fixes)
hist = ds.hist_CHL.squeeze()
hist = hist.sel(zone=ZONE)
# join = ['low' if m == 'low' else 'hi' for m in hist.mask]
# hist = xr.core.groupby.DataArrayGroupBy(
# hist, xr.DataArray(join, dims=['mask'], name='mask')
# ).sum(keep_attrs=True)
def filter(da, N=2, Wn=0.1, **kwargs):
params = signal.butter(N, Wn, **kwargs)
filt = xr.apply_ufunc(
lambda x: signal.filtfilt(*params, x), da,
input_core_dims=[['time']],
output_core_dims=[['time']],
dask='parallelized',
output_dtypes=[('float')]
)
return filt
colors = plot_util.get_mask_colors()
def find_peak(hist):
year = str(int(hist.time.dt.year[0]))
print(year)
hist = hist.sel(time=slice(year+DATE_START, year+DATE_STOP))
if DAY_RES > 1:
hist = hist.resample(time='{}D'.format(DAY_RES)).sum(keep_attrs=True)
hist = hist.where(hist.sum('bins_CHL') > 1500)
pdf = lh.get_pdf(hist)
med = lh.get_median(pdf).load()
# Remove NaNs
med = med.interpolate_na('time', fill_value='extrapolate')
# Sampling freq (in day^-1)
fs = 1./1
# Cutoff frequency
fcut = 1./25
Wn = fcut/(fs/2)
# Butterworth order
Bn = 2
if METHOD[0] == 'filter':
filt = filter(med, Bn, Wn)
elif METHOD[0] == 'mean':
filt = med.rolling(time=METHOD[1]).mean().dropna('time')
elif METHOD[0] == 'median':
filt = med.rolling(time=METHOD[1]).median().dropna('time')
times = filt.time[filt.argmax('time')]
dchl = filter(med, Bn, Wn).diff('time')
dchl = dchl.isel(time=slice(10, -10))
dchl = dchl.resample(time='3H').interpolate('cubic')
times_dchl = dchl.time[dchl.argmax('time')]
cset = tc.tol_cset('vibrant')
m.sel(mask='bkg').plot(ax=ax, color=cset.blue)
m.sel(mask='frt').plot(ax=ax, color=cset.red)
ax.axvline(tmax['bkg'][i], color=cset.blue)
ax.axvline(tmax['frt'][i], color=cset.red)
if PLOT:
fig, axes = plt.subplots(2, 1, figsize=(10, 5), num=year, sharex=True)
fig.subplots_adjust(left=0.05, bottom=0.05, right=0.98, top=0.95)
ax1, ax2 = axes
PLOT = True
DAY_RES = 2
DATE_START = '0201'
DATE_STOP = '0601'
for m, c in colors.items():
ax1.plot(med.time, med.sel(mask=m), color=c, lw=0.5)
ax1.plot(filt.time, filt.sel(mask=m), color=c)
ax1.axvline(times.sel(mask=m).values, color=c)
args = lib.get_args(['region', 'days', 'fixes'])
args['fixes'].setdefault('var', 'CHL')
args['fixes'].setdefault('zone', 'GS3_N')
ax2.plot(dchl.time, dchl.sel(mask=m), color=c)
ax2.axvline(times_dchl.sel(mask=m).values, color=c)
ds = lib.data.hists.get_data(args, remove_dims=['threshold'])
ds = ds.squeeze()
ax1.set_xlim(*med.time[[0, -1]].values)
ds = lib.data.hists.resample_8D_multiple_years(ds, days=DAY_RES)
ds.load()
lib.data.hists.normalize_hist(ds)
med = lib.data.hists.get_percentile(ds, 0.5)
out = xr.concat([t.squeeze().reset_coords(drop=True)
for t in [times, times_dchl]],
dim=xr.DataArray(['max', 'dchl'], dims=['kind']))
return out.squeeze().reset_coords(drop=True)
years = list(set(ds.time.dt.year.values))
tmax = dict(bkg=[], frt=[])
for year in years:
ys = '{:04d}'.format(year)
m = med.sel(time=slice(ys+DATE_START, ys+DATE_STOP))
for mask in m.mask.values:
i = m.sel(mask=mask).argmax()
tmax[mask].append(m.time[i].values)
tmax = hist.groupby('time.year').map(find_peak)
days = tmax.dt.dayofyear
if PLOT:
plot(med.sel(time=slice(ys+'0101', ys+'1231')), tmax)
colors = plot_util.get_mask_colors()
colors.pop('low')
fig, ax = plt.subplots(figsize=(6, 5))
fig.subplots_adjust(left=0.08, bottom=0.05, right=0.98, top=0.98)
fig.subplots_adjust(left=0.12, bottom=0.08, right=0.98, top=0.93)
markers = dict(max='o', dchl='x')
tp = 'datetime64[D]'
tmax_bkg = np.array(tmax['bkg'], tp)
tmax_frt = np.array(tmax['frt'], tp)
delta = tmax_bkg - tmax_frt
for m in ['mid', 'hi']:
for kind, ms in markers.items():
d = days.sel(kind=kind)
ax.scatter(d.sel(mask='low'), d.sel(mask=m),
c=colors[m], marker=ms)
ax.scatter(years, delta)
for year in [2003, 2010, 2014]:
for m in ['mid', 'hi']:
for kind, ms in markers.items():
invalid = plt.matplotlib.markers.MarkerStyle(
ms, fillstyle='none')
d = days.sel(year=year, kind=kind)
ax.scatter(d.sel(mask='low'), d.sel(mask=m),
color='orange', marker=invalid)
pos = years[-1] + 1
ax.boxplot(delta, positions=[pos],
meanline=True, manage_ticks=False, widths=0.7)
# ax.yaxis.set_major_locator(plt.matplotlib.ticker.MultipleLocator(DAY_RES))
ax.set_xlabel('Day of max Chl-a for low')
ax.set_ylabel('Day of max Chl-a for mid, hi')
ax.set_aspect('equal')
ax.set_ylabel('Avance des fronts [Jours]')
c = int(days.mean())
ax.axline((c, c), slope=1, color='k')
fig.canvas.draw()
for tick in ax.xaxis.get_major_ticks():
if (lab := tick.label1).get_text() == str(pos):
lab.set_visible(False)
fig.legend(*plot_util.get_mask_legend(colors),
bbox_to_anchor=[0.98, 0.99], loc='upper right', ncol=len(colors),
borderaxespad=0.)
fig.savefig(path.join(lib.root_plot, 'Hists', 'bloom_shift',
'fig_{}_{}.png'.format(*METHOD)))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment