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

Add HF_inv scripts

dataset for inversion dates
new argument: season
parent c50ee7d2
"""Find and write heat fluxes inversions in boxes."""
from os import path
import matplotlib.pyplot as plt
import numpy as np
from scipy import ndimage
import tol_colors as tc
import lib
import lib.data.HF_TS
import lib.data.hists
# plt.switch_backend('agg')
cset = tc.tol_cset('vibrant')
def filter(da):
return da.rolling(dict(time=8), center=True).median()
def set_ax_limit(da, ax):
vmax = da.max()
vmin = da.min()
margin = (vmax-vmin)*0.1
ax.set_ylim(vmin-margin, vmax+margin)
def find_start_idx(hf, steps_threshold=8):
cond = hf > 0
regions = ndimage.find_objects(ndimage.label(cond)[0])
for slc in regions:
s = np.sum(cond[slc])
if s >= steps_threshold:
return slc[0].start + steps_threshold
raise RuntimeError(f"No {steps_threshold} consecutive positive values")
def add_args(parser):
parser.add_argument('-plot', action='store_true')
args = lib.get_args(['region', 'year', 'fixes', 'grid_file'],
add_args=add_args)
args['fixes']['Y'] = args['year']
PLOT = args['plot']
args['fixes'].update(dict(
number=2,
coef=0,
scale=30,
threshold=15.0,
zone=args['grid_file']
))
hf = lib.data.HF_TS.get_data(args)
if PLOT:
hist = lib.data.hists.get_data(args)
hist = hist.sel(variable='CHL').squeeze().load()
hist_ = hist['hist'].resample(time='8D').sum().to_dataset()
hist_['bins'] = hist['bins']
hist = hist_.load()
hist = lib.data.hists.normalize_hist(hist)
med = lib.data.hists.get_percentile(hist, .5)
med = med.where(hist.n.sel(mask='frt') > 5)
if PLOT:
fig, ax = plt.subplots(figsize=(10, 5))
fig.subplots_adjust(left=0.10, bottom=0.05, right=0.95, top=0.94)
axb = ax.twinx()
hf['hf'] = hf.hf / (3600*24)
ax.axhline(0, color='k')
_hf_ts, = ax.plot(hf.time, hf.hf.isel(idx=0), ls=':')
_hf_f, = ax.plot(hf.time, filter(hf.hf.isel(idx=0)))
_med_bkg, = axb.plot(med.time, med.isel(zone=0).sel(mask='bkg'),
ls='-', color=cset.teal)
_med_frt, = axb.plot(med.time, med.isel(zone=0).sel(mask='frt'),
ls='-', color=cset.orange)
_vline = ax.axvline(hf.time[0].values, color='k', ls='--')
ax.set_ylabel("Heat Surface Flux [W/m²]")
axb.set_ylabel("Chl-a Front/Background")
ax.set_xlim(hf.time[0], hf.time[-1])
inv_times = []
for idx in hf.idx.values:
hfi = hf.sel(idx=idx)
hfi_f = filter(hfi.hf).load()
if PLOT:
medi = med.sel(zone=idx)
ax.set_title(idx)
_hf_ts.set_data(hfi.time, hfi.hf)
_hf_f.set_data(hfi.time, hfi_f)
_med_bkg.set_data(medi.time, medi.sel(mask='bkg'))
_med_frt.set_data(medi.time, medi.sel(mask='frt'))
set_ax_limit(hfi.hf, ax)
set_ax_limit(med, axb)
try:
inv_idx = find_start_idx(hfi_f.to_masked_array().data)
inv_time = str(hfi_f.time[inv_idx].dt.strftime('%F').values)
except RuntimeError:
inv_idx = 0
inv_time = ''
if PLOT:
_vline.set_data(np.repeat(hfi_f.time[inv_idx].values, 2), (0, 1))
inv_times.append((idx, inv_time))
if PLOT:
fig.canvas.draw()
fig.savefig(path.join(lib.root_plot, 'HF', 'TS',
'find_inversion_{}.png'.format(idx)),
dpi=150)
ofile = path.join(lib.root_data, args['region'], 'HF',
args['grid_file'] + '_inversion_times_{d}.txt'.format(
args['year']))
with open(ofile, 'w') as f:
for t in inv_times:
print(*t, file=f)
from datetime import datetime
from os import path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import tol_colors as tc
import xarray as xr
import lib
import lib.data.frt_TS
cset = tc.tol_cset('vibrant')
args = lib.get_args(['region', 'grid_file', 'year', 'fixes',
'scale', 'number', 'coef', 'threshold',
'season'])
# Parse HF inversion dates
inv_file = path.join(lib.root_data, args['region'], 'HF',
args['grid_file'] + '_inversion_times_{}.txt'.format(
args['year']))
inv_times = {}
with open(inv_file, 'r') as infile:
for line in infile.readlines():
idx, date = line.split()
inv_times[idx] = date
# Load Chl TS data
args['zone'] = args['grid_file']
args['fixes']['Y'] = args['year']
args['days'] = 8
chl = lib.data.frt_TS.get_data(args).sel(variable='CHL')
chl = chl.sel(time=slice(datetime(args['year'], 1, 1),
datetime(args['year'], 6, 30)))
dchl = chl.med.diff('time')
dchl = dchl.rename(zone='idx')
"""
idx = '0602'
inv_time = datetime.strptime(inv_times[idx], '%Y-%m-%d')
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
fig.subplots_adjust(left=0.05, bottom=0.10, right=0.98, top=0.98)
ax1, ax2 = axes
dchl.sel(idx=idx, mask='bkg').plot.line(ax=ax1, color=cset.teal)
dchl.sel(idx=idx, mask='frt').plot.line(ax=ax1, color=cset.orange)
ax1.axvline(inv_time, color='r', ls='--')
lag = pd.to_datetime(dchl.time.values) - inv_time
dchl_lag = dchl.assign_coords(time=lag.days)
dchl_lag = dchl_lag.rename(time='lag')
dchl_lag.lag['units'] = 'days since {}'.format(inv_time.strftime('%F'))
dchl_lag.sel(idx=idx, mask='bkg').plot.line(ax=ax2, color=cset.teal)
dchl_lag.sel(idx=idx, mask='frt').plot.line(ax=ax2, color=cset.orange)
"""
idx_sel = ['0601', '0602', '0701', '0702']
dchl_lag_sel = []
for idx in idx_sel:
inv_time = datetime.strptime(inv_times[idx], '%Y-%m-%d')
lag = pd.to_datetime(dchl.time.values) - inv_time
dchl_lag = dchl.sel(idx=[idx]).assign_coords(time=lag.days)
dchl_lag = dchl_lag.rename(time='lag')
dchl_lag.lag['units'] = 'days'
dchl_lag_sel.append(dchl_lag)
s = xr.concat(dchl_lag_sel, dim='idx', join='outer', fill_value=np.nan)
fig, ax = plt.subplots(figsize=(6, 5))
fig.subplots_adjust(left=0.15, bottom=0.10, right=0.98, top=0.98)
for ds in dchl_lag_sel:
ax.plot(ds.lag, ds.sel(mask='bkg'), lw=0.7)
ax.plot(s.lag, s.sel(mask='bkg').mean('idx'), color='k')
ax.set_xlabel('Jours depuis inversion flux chaleur')
ax.set_ylabel('dChl (background)')
......@@ -123,6 +123,7 @@ def get_args(args, description='', add_args=None):
add_arg('zone', str)
add_arg('level', str)
add_arg('grid_file', str)
add_arg('season', str)
if 'fixes' in args:
parser.add_argument('-fix', nargs=2, type=str, action='append',
......
......@@ -10,3 +10,4 @@ coef: 0
level: L3
zone: gyre
grid_file: boxgrid_5.0_5.0
season: spring
"""Date of Surface Heat Fluxes inversion."""
from datetime import datetime
from os import path
import lib
import lib.data
ARGS_DIR = {'region', 'season'}
defaults = {'season': 'spring'}
def get_root(args=None, **kwargs):
args = lib.data.process_args(ARGS_DIR, args, replace_defaults=defaults,
**kwargs,)
root = path.join(lib.root_data, args['region'], 'HF', 'inversion',
args['season'])
return root
auto_attr = lib.data.create_data(__name__, pregex, get_root, ARGS_DIR,
defaults)
def get_data(args=None, to_datetime=False, **kwargs):
finder = get_finder(args, **kwargs)
inv_times = {}
for file in finder.get_files():
with open(file, 'r') as f:
for line in f.readlines():
idx, date = line.split()
year = int(date[:4])
if to_datetime:
date = datetime.strptime(date, '%Y-%m-%d')
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