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

Up plots

parent 5cac9485
......@@ -2,7 +2,7 @@
from os import path
import matplotlib.pyplot as plt
import tol_colors as tc
import xarray as xr
import lib
import lib.data.hists as lh
......@@ -10,36 +10,51 @@ import Plots.util as plot_util
plot_util.use_tex(True)
args = dict(region='GS', days=1, kind='2thr', mask=True)
args = dict(region='GS', days=1, kind='2thr')
fixes = dict(
scale=30.,
number=2,
coef=0,
thr_lo=5.,
thr_hi=15.
thr_hi=10.
)
ds = lh.get_data(args, fixes=fixes)
ds = lh.get_data(args, fixes=fixes, mask=True)
ds = ds.squeeze()
ds = lh.resample_8D_multiple_years(ds)
ds = ds.load()
# join = ['hi' if m == 'low' else 'low' for m in ds.mask]
# ds = xr.core.groupby.DatasetGroupBy(
# ds, xr.DataArray(join, dims=['mask'], name='mask')
# ).sum(keep_attrs=True)
pdf = lh.get_pdf(ds)
ds['med'] = lh.get_median(pdf).load()
N = lh.get_N(ds.hist_CHL)
ds = ds.drop_vars(['hist_CHL', 'hist_sst'])
ds_nm = lh.get_data(args, fixes=fixes, mask=False)
ds_nm = ds_nm.squeeze()
ds_nm = lh.resample_8D_multiple_years(ds_nm)
ds_nm = ds_nm.load()
N = lh.get_N(ds_nm.hist_CHL)
frac = N / N.sum('mask') * 100.
ds['frac'] = frac
pdf = lh.get_pdf(ds)
ds['med'] = lh.get_median(pdf).load()
T = lh.get_total(ds_nm.hist_CHL)
ds = ds.drop_vars(['hist_CHL', 'hist_sst'])
def get_C_mean(T, N, masks):
Ts, Ns = [da.sel(mask=masks).sum('mask') for da in [T, N]]
return Ts/Ns
# Chl moyenne
C = get_C_mean(T, N, ['low', 'mid', 'hi'])
C_nf_hi = get_C_mean(T, N, ['low', 'mid'])
C_nf_mid = get_C_mean(T, N, ['low', 'hi'])
C_nf = xr.concat([C_nf_mid, C_nf_hi], dim=xr.DataArray(['mid', 'hi'],
dims=['mask']))
ds['surplus'] = (C-C_nf)/C_nf * 100.
clis = lh.get_climatology(ds)
clis = lh.get_climatology(ds.load())
def loop_zones():
......@@ -57,10 +72,10 @@ def plot(cli, cli_d1, cli_d9, ax, var, sel, color='k', ls='-'):
color=color, alpha=.3)
colors = plot_util.get_mask_colors('vibrant', 'red', 'blue', 'teal')
colors = plot_util.get_mask_colors()
fig, axes = plt.subplots(3, 3, figsize=(10, 6), sharex=True)
fig.subplots_adjust(left=0.06, bottom=0.05, right=0.98, top=0.95,
fig, axes = plt.subplots(4, 3, figsize=(8, 6), sharex=True)
fig.subplots_adjust(left=0.08, bottom=0.05, right=0.98, top=0.95,
hspace=0.1, wspace=0.15)
for ax in axes.flat:
......@@ -80,8 +95,10 @@ for j, zone in enumerate(loop_zones()):
for m in ['mid', 'hi']:
sel['mask'] = m
plot(*clis, axes[2, j], 'frac', sel, color=colors[m])
plot(*clis, axes[3, j], 'surplus', sel, color=colors[m])
axes[2, j].set_ylim(0, 60)
plt.setp(axes[2, j].xaxis.get_majorticklabels(), ha='left')
axes[3, j].set_ylim(-6, 23)
plt.setp(axes[3, j].xaxis.get_majorticklabels(), ha='left')
......@@ -90,6 +107,7 @@ axes[0, 0].set_xlim(*(clis[0].time[[0, -1]].values))
axes[0, 0].set_ylabel(plot_util.SST_LABEL)
axes[1, 0].set_ylabel(plot_util.CHL_LABEL)
axes[2, 0].set_ylabel('Fronts surface [\%]')
axes[3, 0].set_ylabel('Surplus [\%]')
axes[0, 0].set_title('South')
axes[0, 1].set_title('Jet')
......
"""Sensitivity to parameters.
4 subplots. Each showing Chl histogram for a full year in background and front
for different values of parameters. One parameter is changed each time. The
threshold for each parameter is choosen so that the max surface fraction of
fronts is ~10%.
2 tests: varying rolling window size (scale) and coefficients weights.
Each test is done for all 3 GS3 zones.
Take zones south of jet (GS3_I and GS3_S).
Threshold is not changed between tests.
Test only done with ostia data.
"""
from os import path
......@@ -38,34 +37,13 @@ fixes_def = dict(
def get_data(fixes):
def preprocess(ds):
d = 'threshold'
if d in ds.coords and ds.coords[d].size <= 1:
ds = ds.isel({d: 0})
return ds.reset_coords([d], drop=True)
ds = lh.get_data(args, fixes=fixes, open_mf_kw=dict(process=preprocess))
if 'CHL' in ds.VARS:
var = 'CHL'
else:
var = 'chl_ocx'
hist = ds[lh.var_name('hist', var)].squeeze()
# if 'zone' in hist.dims:
# hist = hist.sel(zone=fixes['zone'])
ds = lh.get_data(args, fixes=fixes)
hist = ds.hist_CHL.squeeze()
hist = hist.sum('time', keep_attrs=True)
hist.load()
pdf = lh.get_pdf(hist)
med = lh.get_median(pdf)
return med, pdf, pdf[lh.var_name('bins', var)]
def loop_thr(fixes_iter):
for values in zip(*fixes_iter.values()):
fixes = fixes_def.copy()
for k, v in zip(fixes_iter.keys(), values):
fixes[k] = v
yield fixes
return med, pdf, pdf.bins_CHL
cset = tc.tol_cset('high-contrast')
......@@ -81,71 +59,59 @@ med_bkg_tot = []
## Subplot Scale
print("Scale", flush=True)
colors = [cset.yellow, cset.red, cset.blue]
med_frt = []
med_bkg = []
fixes_iter = dict(
scale=[20., 30., 40.],
threshold=[5., 5., 5.]
)
fixes = fixes_def.copy()
fixes['scale'] = [20., 30., 40.]
med, pdf, bins = get_data(fixes)
labels = []
handles = []
for i, fixes in enumerate(loop_thr(fixes_iter)):
med, pdf, bins = get_data(fixes)
for z, ax in zip('NIS', axes[:, 0]):
zone = 'GS3_' + z
ax.plot(bins, pdf.sel(mask='frt', zone=zone), color=colors[i], ls='-')
ax.plot(bins, pdf.sel(mask='bkg', zone=zone), color=colors[i], ls='--')
med_frt.append(med.sel(mask='frt', zone=zone).values)
med_bkg.append(med.sel(mask='bkg', zone=zone).values)
ax.axvline(med_frt[-1], color=colors[i], ls='-')
ax.axvline(med_bkg[-1], color=colors[i], ls='--')
for z, ax in zip('NIS', axes[:, 0]):
zone = 'GS3_' + z
for i, s in enumerate(fixes['scale']):
sel = dict(zone=zone, scale=s)
ax.plot(bins, pdf.sel(mask='frt', **sel), color=colors[i], ls='-')
ax.plot(bins, pdf.sel(mask='bkg', **sel), color=colors[i], ls='--')
ax.axvline(med.sel(mask='frt', **sel), color=colors[i], ls='-')
ax.axvline(med.sel(mask='bkg', **sel), color=colors[i], ls='--')
labels.append(r'\SI{{{:.0f}}}{{\kilo\meter}}'.format(fixes['scale']))
handles.append(plt.Line2D((), (), color=colors[i], ls='-'))
labels = [r'\SI{{{:.0f}}}{{\kilo\meter}}'.format(s) for s in fixes['scale']]
handles = [plt.Line2D((), (), color=c) for c in colors]
axes[1, 0].legend(labels=labels, handles=handles, framealpha=1)
print("FRT std: ", np.std(med_frt))
med_frt_tot += list(np.array(med_frt) - np.mean(med_frt))
med_bkg_tot += list(np.array(med_bkg) - np.mean(med_bkg))
med_norm = med - med.mean('scale')
print("FRT std: ", med_norm.sel(mask='frt').std().values)
med_frt_tot += list(med_norm.sel(mask='frt').values.reshape(-1))
med_bkg_tot += list(med_norm.sel(mask='bkg').values.reshape(-1))
## Subplot Variant
print("Variant")
ax = axes[0, 1]
colors = [cset.black, cset.yellow, cset.red, cset.blue]
med_frt = []
med_bkg = []
fixes_iter = dict(
coef=[0, 1, 2, 3],
threshold=[5., 5., 5., 5.]
)
for i, fixes in enumerate(loop_thr(fixes_iter)):
med, pdf, bins = get_data(fixes)
fixes = fixes_def.copy()
fixes['coef'] = [0, 1, 2, 3]
med, pdf, bins = get_data(fixes)
for z, ax in zip('NIS', axes[:, 1]):
zone = 'GS3_' + z
ax.plot(bins, pdf.sel(mask='frt', zone=zone), color=colors[i], ls='-')
ax.plot(bins, pdf.sel(mask='bkg', zone=zone), color=colors[i], ls='--')
med_frt.append(med.sel(mask='frt', zone=zone).values)
med_bkg.append(med.sel(mask='bkg', zone=zone).values)
ax.axvline(med_frt[-1], color=colors[i], ls='-')
ax.axvline(med_bkg[-1], color=colors[i], ls='--')
for z, ax in zip('NIS', axes[:, 1]):
zone = 'GS3_' + z
for i, coef in enumerate(fixes['coef']):
sel = dict(zone=zone, coef=coef)
ax.plot(bins, pdf.sel(mask='frt', **sel), color=colors[i], ls='-')
ax.plot(bins, pdf.sel(mask='bkg', **sel), color=colors[i], ls='--')
ax.axvline(med.sel(mask='frt', **sel), color=colors[i], ls='-')
ax.axvline(med.sel(mask='bkg', **sel), color=colors[i], ls='--')
labels = ['Balanced', 'More skewness', 'More variance', 'More bimodality']
handles = [plt.Line2D((), (), color=c) for c in colors]
axes[1, 1].legend(labels=labels, handles=handles, framealpha=1)
print("FRT std: ", np.std(med_frt))
med_frt_tot += list(np.array(med_frt) - np.mean(med_frt))
med_bkg_tot += list(np.array(med_bkg) - np.mean(med_bkg))
med_norm = med - med.mean('coef')
print("FRT std: ", med_norm.sel(mask='frt').std().values)
med_frt_tot += list(med_norm.sel(mask='frt').values.reshape(-1))
med_bkg_tot += list(med_norm.sel(mask='bkg').values.reshape(-1))
print("TOTAL FRT std: ", np.std(med_frt_tot))
print("TOTAL BKG std: ", np.std(med_bkg_tot))
......
......@@ -25,7 +25,7 @@ fixes = dict(
ds = lh.get_data(args, fixes=fixes)
ds = ds.squeeze()
ds = lh.resample_8D_multiple_years(ds, days=3)
pdf = lh.get_pdf(ds, variables='CHL')
pdf = lh.get_pdf(ds, variables='CHL').load()
med = lh.get_median(pdf)
med = med.load()
......@@ -38,11 +38,23 @@ right = 0.96
fig.subplots_adjust(left=0.08, bottom=0.05, right=right, top=0.93,
hspace=0.1)
med_frt_std = 0.0041
med_bkg_std = 0.0020
error = dict(
low=med_bkg_std*2,
mid=med_frt_std*2,
hi=med_frt_std*2
)
for ax, z in zip(axes, 'NIS'):
zone = lh.get_GS3_name(z)
for m, c in colors.items():
ax.plot(med.time, med.sel(zone=zone, mask=m),
color=c)
M = med.sel(zone=zone, mask=m)
ax.plot(med.time, M, color=c)
ax.fill_between(med.time.values, M-error[m], M+error[m],
color=c, alpha=0.3)
if z != 'N':
ax.set_ylim(None, ymax.sel(zone=zone))
......@@ -52,8 +64,8 @@ for ax, z in zip(axes, 'NIS'):
for i, ax in enumerate(axes):
ax.annotate(chr(97+i)+')', (0., 1.), xycoords='axes fraction',
xytext=(-40, 2), textcoords='offset points',
ha='left', va='bottom')
xytext=(-5, -5), textcoords='offset points',
ha='right', va='bottom')
ax.xaxis.axis_date()
ax.set_ylabel(plot_util.CHL_LABEL, labelpad=10.0)
......
......@@ -16,19 +16,13 @@ plot_util.set_style()
plot_util.use_tex(True)
args = dict(
region='GS',
days=1
)
args = dict(region='GS', days=1, kind='2thr')
zone = 'box7'
fixes = dict(
scale=30.,
number=2,
coef=0,
zone=zone,
threshold=15.0
# threshold=dict(S=6.0, N=15.0)[zone[-1]]
thr_lo=15.0
)
ds = lib.data.hists.get_data(args, fixes=fixes)
......@@ -83,23 +77,20 @@ def plot(ax, var, sel, color='k', ls='-'):
## Plot
fig, axes = plt.subplots(3, 1, figsize=(6, 4), sharex=True)
fig, axes = plt.subplots(4, 1, figsize=(6, 4), sharex=True)
fig.subplots_adjust(left=0.11, bottom=0.07, right=0.98, top=0.935,
hspace=0.1)
ax1, ax2, ax3 = axes
ax1, ax2, ax3, ax4 = axes
cset = tc.tol_cset('vibrant')
plot(ax1, 'med', dict(variable='CHL', mask='frt'),
color=cset.orange)
plot(ax1, 'med', dict(variable='CHL', mask='bkg'),
color=cset.teal)
plot(ax2, 'med', dict(variable='sst', mask='frt'),
color=cset.orange)
plot(ax2, 'med', dict(variable='sst', mask='bkg'),
color=cset.teal)
plot(ax1, 'med', dict(variable='CHL', mask='frt'), color=cset.orange)
plot(ax1, 'med', dict(variable='CHL', mask='bkg'), color=cset.teal)
plot(ax2, 'med', dict(variable='sst', mask='frt'), color=cset.orange)
plot(ax2, 'med', dict(variable='sst', mask='bkg'), color=cset.teal)
plot(ax3, 'frac', dict())
plot(ax4, 'surplus', dict())
labelpad = 11
......
......@@ -15,12 +15,7 @@ plot_util.set_style()
plot_util.use_tex(True)
args = dict(
region='GS',
days=1
)
zone = 'GS3_N'
args = dict(region='GS', days=1, kind='2thr')
fixes = dict(
scale=30.,
number=2,
......@@ -29,11 +24,10 @@ fixes = dict(
var='CHL'
)
start = 2005
# end = start + 4
end = 2015
start = 2007
end = start + 4
ds = lib.data.hists.get_data(args, remove_dims=['threshold'], fixes=fixes)
ds = lib.data.hists.get_data(args, fixes=fixes)
ds = ds.squeeze()
ds = ds.sel(time=slice(datetime(start, 1, 1), datetime(end, 12, 31)))
ds = ds.resample(time='8D').sum()
......
......@@ -7,9 +7,14 @@ import tol_colors as tc
import lib
import lib.data.hists as lh
import Plots.util as plot_util
plot_util.set_style()
plot_util.use_tex(True)
args = lib.get_args(['region', 'days', 'fixes'])
args['kind'] = '2thr'
fixes = dict(
scale=30.,
number=2,
......@@ -22,39 +27,41 @@ args['fixes'] = fixes
ds = lh.get_data(args)
ds = ds.squeeze()
ds = ds.groupby('time.season').sum(keep_attrs=True)
ds.load()
pdf = lh.get_pdf(ds)
cset = tc.tol_cset('vibrant')
colors = plot_util.get_mask_colors()
for var in pdf.VARS:
fig, axes = plt.subplots(3, 4, figsize=(10, 8), sharex=True)
fig.subplots_adjust(left=0.05, bottom=0.05, right=0.98, top=0.94,
fig, axes = plt.subplots(3, 4, figsize=(8, 6), sharex=True)
fig.subplots_adjust(left=0.05, bottom=0.08, right=0.98, top=0.95,
wspace=0.1, hspace=0.1)
fig.suptitle(var)
for i, (z, z_title) in enumerate(zip('NIS', ['North', 'Jet', 'South'])):
axes[i, 0].set_ylabel(z_title)
zone = 'GS3_{}'.format(z)
for j, season in enumerate(['DJF', 'MAM', 'JJA', 'SON']):
axes[0, j].set_title(season)
p = pdf[lh.var_name('pdf', var)].sel(zone=zone, season=season)
for mask, color in zip(['low', 'mid', 'hi'],
['cyan', 'orange', 'red']):
for mask, color in colors.items():
axes[i, j].plot(p[lh.var_name('bins', var)],
p.sel(mask=mask), color=getattr(cset, color))
p.sel(mask=mask), color=color)
for ax in axes.flat:
ax.tick_params(axis='y', labelleft=False)
for ax in axes[-1, :]:
ax.set_xlabel('--')
if var == 'CHL':
axes[-1, 0].set_xlabel(plot_util.CHL_LABEL)
axes[0, 0].set_xlim(0, 1.1)
elif var == 'sst':
axes[-1, 0].set_xlabel(plot_util.SST_LABEL)
axes[0, 0].set_xlim(-4, 32)
axes[1, 2].legend(*plot_util.get_mask_legend(colors))
fig.canvas.draw()
fig.savefig(path.join(lib.root_plot, 'Hists', 'hists_season',
'fig_{}.png'.format(var)),
......
......@@ -14,7 +14,6 @@ fixes = dict(
scale=30.,
number=2,
coef=0,
zone='GS3_N',
var='CHL'
)
......
......@@ -61,7 +61,7 @@ fig.delaxes(axes[0, -1])
def add_cbar(im, ax):
bbox = ax.bbox.transformed(fig.transFigure.inverted())
cax = fig.add_axes([bbox.xmax, bbox.ymin, 0.01, bbox.ymax-bbox.ymin])
cax = fig.add_axes([bbox.xmax, bbox.ymin, 0.01, bbox.height])
cbar = fig.colorbar(im, cax=cax)
return cbar
......
import itertools
from os import path
import cartopy.crs as ccrs
......@@ -10,10 +11,12 @@ from matplotlib.colors import LogNorm
import lib
import lib.data.ostia as lo
import lib.data.globcolour as lg
import lib.data.SN_separation as lthr
import Plots.util as plot_util
import Compute.compute_temp_separation_zones as csep
plot_util.set_style()
plot_util.use_tex(True)
args = dict(
region='GS',
level='L3',
......@@ -36,33 +39,98 @@ def normalize_hist(bins, hist):
sst = lo.get_data(args).sst.squeeze()
chl = lg.get_data(args).CHL.squeeze()
sep = csep._get_separation(sst.sel(lat=slice(None, 32)).to_numpy())
thr = sep['threshold']
## Plot
var_range = dict(sst=[-4, 30], CHL=[1e-1, 5])
var_range_kw = {v: dict(zip(['vmin', 'vmax'], r)) for v, r in var_range.items()}
cset_sep = tc.tol_cset('bright')
zone_names = dict(N='North', I='Jet', S='South')
fig = plt.figure(figsize=(7.6, 6))
gs = fig.add_gridspec(2, 1, left=0.04, bottom=0.07, right=0.94, top=0.97,
hspace=0.1)
## Maps
axes_map = gs[0, 0].subgridspec(1, 2, wspace=0.2).subplots(
subplot_kw=dict(projection=ccrs.PlateCarree()))
axes_map[0].sharex(axes_map[1])
axes_map[0].sharey(axes_map[1])
mult_sst = 9./0.150 # So that hist is 9 at max
im_sst = sst.plot.imshow(
ax=axes_map[0], cmap=cm.thermal, **var_range_kw['sst'],
add_labels=False, add_colorbar=False
)
im_chl = chl.plot.imshow(
ax=axes_map[1], cmap=cm.algae, **var_range_kw['CHL'],
add_labels=False, add_colorbar=False, norm=LogNorm()
)
for ax in axes_map:
ax.xaxis.set_major_formatter(plot_util.LonFormatter('.0f'))
ax.yaxis.set_major_formatter(plot_util.LatFormatter('.0f'))
ax.set_aspect('equal')
fig, axes = plt.subplots(2, 2, figsize=(8, 6))
fig.subplots_adjust(left=0.06, bottom=0.05, right=0.97, top=0.95)
ax.coastlines(color='k', lw=0.7)
gl = ax.gridlines(draw_labels=True, lw=0.7)
gl.top_labels = False
gl.right_labels = False
gl.xformatter = plot_util.LonFormatter('.0f')
gl.yformatter = plot_util.LatFormatter('.0f')
sst.plot.contour(ax=ax, levels=[thr], colors=cset_sep.purple,
add_labels=False)
ax.axhline(32, color=cset_sep.purple, ls='--')
axes[0, 0].sharex(axes[0, 1])
axes[0, 0].sharey(axes[0, 1])
sst.plot.imshow(ax=axes[0, 0], cmap=cm.thermal, **var_range_kw['sst'],
cbar_kwargs=dict(label='SST', pad=0.02))
chl.plot.imshow(ax=axes[0, 1], cmap=cm.algae, **var_range_kw['CHL'],
cbar_kwargs=dict(label='CHL', pad=0.02), norm=LogNorm())
ax = axes_map[0]
kw = dict(color='#E2E2E2', ha='center', size='large')
lat_str = r'\textbf{{{}}}'
ax.annotate(lat_str.format(zone_names['N']), (-46, 48), **kw)
ax.annotate(lat_str.format(zone_names['I']), (-50, 35), **kw)
ax.annotate(lat_str.format(zone_names['S']), (-55, 24), **kw)
var_labels = [plot_util.SST_LABEL, plot_util.CHL_LABEL]
caxes = []
for ax, im, lab in zip(axes_map, [im_sst, im_chl], var_labels):
cax = plot_util.add_cax(ax, axes_class=plt.Axes)
plt.colorbar(im, cax=cax, label=lab)
caxes.append(cax)
## Hists
fig.canvas.draw()
# Impose positions of axes from axes_map
axes_hist = []
for ax, cax in zip(axes_map, caxes):
bbox = ax.bbox.transformed(fig.transFigure.inverted())
bbox_cax = cax.bbox.transformed(fig.transFigure.inverted())
axes_hist.append(fig.add_axes(
(bbox.xmin, gs.bottom,