Skip to content

GitLab

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

Up plots

parent eba8c7d1
Hide whitespace changes
Inline Side-by-side
Showing with 288 additions and 460 deletions
Plots/Hists/all_zones_median_cli.py
View file @ 5cac9485
......@@ -10,13 +10,13 @@ import Plots.util as plot_util
plot_util.use_tex(True)
args = dict(region='GS', days=1, kind='2thr', mask=False)
args = dict(region='GS', days=1, kind='2thr', mask=True)
fixes = dict(
scale=30.,
number=2,
coef=0,
thr_lo=5.,
thr_hi=10.
thr_hi=15.
)
......@@ -30,9 +30,10 @@ ds = ds.load()
# ds, xr.DataArray(join, dims=['mask'], name='mask')
# ).sum(keep_attrs=True)
N = lh.get_N(ds).sel(var='CHL').reset_coords(drop=True)
N = lh.get_N(ds.hist_CHL)
frac = N / N.sum('mask') * 100.
ds['frac'] = frac
pdf = lh.get_pdf(ds)
ds['med'] = lh.get_median(pdf).load()
......
Plots/Hists/article_hist_sensitivity.py
View file @ 5cac9485
......@@ -17,7 +17,7 @@ import numpy as np
import tol_colors as tc
import Plots.util as plot_util
import lib.data.hists
import lib.data.hists as lh
plot_util.set_style()
plot_util.use_tex(True)
......@@ -25,197 +25,152 @@ plot_util.use_tex(True)
odir = path.join(plot_util.root_plot, 'article_hist_sensitivity')
args = dict(
region='GS',
days=1,
)
args = dict(region='GS', days=1, mask=False, kind='1thr')
fixes = dict(
zone='GS3_(I|S)',
# Default fixes
fixes_def = dict(
Y=2007,
threshold=5.,
number=2,
var='(CHL|chl_ocx)',
Y=2007
scale=30.,
coef=0,
)
sel = dict()
def get_data(fixes, remove_threshold=False):
remove_dims = []
if remove_threshold:
remove_dims.append('threshold')
ds = lib.data.hists.get_data(args, remove_dims=remove_dims, fixes=fixes)
if 'CHL' in ds['variable']:
ds = ds.sel(variable='CHL')
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:
ds = ds.sel(variable='chl_ocx')
ds = ds.squeeze()
ds = ds.sum('time')
ds.load()
# ds['hist'] = ds.hist.rolling(dict(nbins=15), min_periods=1).mean()
ds = lib.data.hists.normalize_hist(ds, True)
med = lib.data.hists.get_percentile(ds, 0.5).persist()
hist = ds['hist'] * 30
return med, hist, hist.bins
var = 'chl_ocx'
hist = ds[lh.var_name('hist', var)].squeeze()
# if 'zone' in hist.dims:
# hist = hist.sel(zone=fixes['zone'])
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)]
cset = tc.tol_cset('high-contrast')
colors = [cset.yellow, cset.red, cset.blue, cset.black]
fig, axes = plt.subplots(2, 2, figsize=(6, 5), sharex=True, sharey=True)
fig.subplots_adjust(left=0.07, bottom=0.11, right=0.97, top=0.90,
wspace=0.15)
med_frt = []
med_bkg = []
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
## Subplot Threshold
print("Threshold")
ax = axes[0, 0]
colors = [cset.yellow, cset.red, cset.blue]
if fixes['zone'].endswith('N'):
thresholds = [10., 12., 15.]
else:
thresholds = [4., 5., 6.] # South
cset = tc.tol_cset('high-contrast')
fixes['coef'] = 0
fixes['scale'] = 30.
fig, axes = plt.subplots(3, 2, figsize=(7, 6))
fig.subplots_adjust(left=0.05, bottom=0.08, right=0.96, top=0.90,
wspace=0.15, hspace=0.15)
labels = []
handles = []
for i, thr in enumerate(thresholds):
fixes['threshold'] = thr
med, hist, bins = get_data(fixes)
ax.plot(bins, hist.sel(mask='frt'), color=colors[i], ls='-')
ax.plot(bins, hist.sel(mask='bkg'), color=colors[i], ls='--')
med_frt.append(med.sel(mask='frt').values)
med_bkg.append(med.sel(mask='bkg').values)
ax.axvline(med_frt[-1], color=colors[i], ls='-')
ax.axvline(med_bkg[-1], color=colors[i], ls='--')
labels.append('threshold: {:.1f}'.format(thr))
handles.append(plt.Line2D((), (), color=colors[i], ls='-'))
med_frt_tot = []
med_bkg_tot = []
ax.legend(labels=labels, handles=handles, framealpha=1)
## Subplot Scale
print("Scale", flush=True)
ax = axes[1, 0]
colors = [cset.yellow, cset.red, cset.blue]
med_frt = []
med_bkg = []
if fixes['zone'] == 'GS_S':
thresholds = [6.0, 7.5, 8.7] # South
else:
thresholds = [15., 17., 19.] # North
scales = [30., 40., 50.]
fixes['coef'] = 0
fixes_iter = dict(
scale=[20., 30., 40.],
threshold=[5., 5., 5.]
)
labels = []
handles = []
for i, (scale, thr) in enumerate(zip(scales, thresholds)):
fixes['scale'] = scale
fixes['threshold'] = thr
med, hist, bins = get_data(fixes)
ax.plot(bins, hist.sel(mask='frt'), color=colors[i], ls='-')
ax.plot(bins, hist.sel(mask='bkg'), color=colors[i], ls='--')
med_frt.append(med.sel(mask='frt').values)
med_bkg.append(med.sel(mask='bkg').values)
ax.axvline(med_frt[-1], color=colors[i], ls='-')
ax.axvline(med_bkg[-1], color=colors[i], ls='--')
labels.append(r'\SI{{{:.0f}}}{{\kilo\meter}}'.format(scale))
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='--')
labels.append(r'\SI{{{:.0f}}}{{\kilo\meter}}'.format(fixes['scale']))
handles.append(plt.Line2D((), (), color=colors[i], ls='-'))
ax.legend(labels=labels, handles=handles, framealpha=1)
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))
## Subplot Variant
print("Variant")
ax = axes[0, 1]
colors = [cset.black, cset.yellow, cset.red, cset.blue]
med_frt = []
med_bkg = []
coefs = list(range(4))
scale = 30.
if fixes['zone'] == 'GS_S':
thresholds = [6.0, 10., 6.0, 10.] # South
else:
thresholds = [15.0, 25., 17.0, 22.] # North
fixes['scale'] = scale
labels = []
handles = []
for i, (coef, thr) in enumerate(zip(coefs, thresholds)):
fixes['coef'] = coef
fixes['threshold'] = thr
med, hist, bins = get_data(fixes)
ax.plot(bins, hist.sel(mask='frt'), color=colors[i], ls='-')
ax.plot(bins, hist.sel(mask='bkg'), color=colors[i], ls='--')
med_frt.append(med.sel(mask='frt').values)
med_bkg.append(med.sel(mask='bkg').values)
ax.axvline(med_frt[-1], color=colors[i], ls='-')
ax.axvline(med_bkg[-1], color=colors[i], ls='--')
if i == 0:
labels.append(r'Base coefficients')
else:
labels.append(r'Variation {:d}'.format(i))
handles.append(plt.Line2D((), (), color=colors[i], ls='-'))
ax.legend(labels=labels, handles=handles, framealpha=1)
fixes_iter = dict(
coef=[0, 1, 2, 3],
threshold=[5., 5., 5., 5.]
)
## Subplot Variant
print("Resolution")
ax = axes[1, 1]
colors = [cset.yellow, cset.red]
for i, fixes in enumerate(loop_thr(fixes_iter)):
med, pdf, bins = get_data(fixes)
fixes['coef'] = 0
fixes['scale'] = 30.
numbers = [1, 2]
thresholds = [2.6, 6.0]
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='--')
labels = []
handles = []
for i, (number, thr) in enumerate(zip(numbers, thresholds)):
fixes['number'] = number
fixes['threshold'] = thr
med, hist, bins = get_data(fixes)
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='--')
ax.plot(bins, hist.sel(mask='frt'), color=colors[i], ls='-')
ax.plot(bins, hist.sel(mask='bkg'), 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)
med_frt.append(med.sel(mask='frt').values)
med_bkg.append(med.sel(mask='bkg').values)
ax.axvline(med_frt[-1], color=colors[i], ls='-')
ax.axvline(med_bkg[-1], color=colors[i], ls='--')
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))
labels.append(r'Resolution \si{{{:d}}}{{km}}'.format([1, 4][i]))
handles.append(plt.Line2D((), (), color=colors[i], ls='-'))
ax.legend(labels=labels, handles=handles, framealpha=1)
print("TOTAL FRT std: ", np.std(med_frt_tot))
print("TOTAL BKG std: ", np.std(med_bkg_tot))
print("BKG std: ", np.std(med_bkg))
print("FRT std: ", np.std(med_frt))
## Cosmetics
for ax in axes[0, :]: # North
ax.set_xlim(0., 0.8)
for ax in axes[1, :]: # Jet
ax.set_xlim(0., 0.5)
for ax in axes[2, :]: # South
ax.set_xlim(0., 0.25)
## Cosmetics
for i, ax in enumerate(axes.flat):
if fixes['zone'] == 'GS_S':
ax.set_xlim(0., 0.25) # South
else:
ax.set_xlim(0., 1.4) # North
titles = ['Scale', 'Coefficients']
for ax, title in zip(axes[0], titles):
ax.set_title(title)
ax.annotate(chr(97+i)+')', (0, 1.0), xycoords='axes fraction',
xytext=(-5, 5), textcoords='offset points',
ha='right', va='bottom')
zones = ['North', 'Jet', 'South']
for ax, z in zip(axes[:, 1], zones):
ax.annotate(z, (1., 0.5), xycoords='axes fraction',
xytext=(5, 0), textcoords='offset points',
ha='left', va='center', rotation='vertical')
for ax in axes[1]:
ax.set_xlabel(r'Chl-\emph{a} [\si{\milli\gram\per\cubic\meter}]')
for ax in axes[-1]:
ax.set_xlabel(plot_util.CHL_LABEL)
fig.canvas.draw()
......@@ -223,13 +178,13 @@ labels = ['Front', 'Background']
handles = [plt.Line2D((), (), color='k', ls=ls)
for ls in ['-', '--']]
bbox = axes[0, 0].bbox.transformed(fig.transFigure.inverted())
fig.legend(labels=labels, handles=handles,
bbox_to_anchor=(bbox.xmin, bbox.ymax+0.01,
bbox.width,
1.-fig.subplotpars.top+0.02),
bbox_to_anchor=(fig.subplotpars.left,
fig.subplotpars.top,
fig.subplotpars.right-fig.subplotpars.left,
1.-fig.subplotpars.top),
loc='lower center', frameon=False,
ncol=2, borderaxespad=0)
ncol=1, borderaxespad=0)
## Savefig
......
Plots/Images/article_example_zoom.py
View file @ 5cac9485
from os import path
from cmocean import cm
from matplotlib import pyplot as plt
from matplotlib.ticker import MultipleLocator
import matplotlib.pyplot as plt
import numpy as np
from cmocean import cm
import lib
import lib.data.ostia
import lib.data.globcolour
import lib.data.modis
import lib.data.hi
import Plots.util as plot_util
from Plots.Images.examples import examples
from lib.box import Box
plot_util.set_style()
plot_util.use_tex(True)
MODIS = False
THRESHOLD = 10.
args = dict(
region='GS',
days=1,
scale=30.,
number=2-int(MODIS),
coef=0
)
ex_id = 1
ex = examples[ex_id]
# Box to show scale
win_frac_bl = [0.3, 0.3]
win_bl = [ex.corner_bl[i] + (ex.corner_tr[i]-ex.corner_bl[i])
* win_frac_bl[i] for i in range(2)]
km_per_deg = 6.371*1e3 * np.pi / 180
win_tr = [bl + args['scale'] / km_per_deg for bl in win_bl]
win = Box.FromCorners(*win_bl, *win_tr)
import lib.data.ostia as lo
import lib.data.hi as lhi
from Plots.Images.examples import Example
fixes = dict(Y=ex.date.year, m=ex.date.month, d=ex.date.day)
THRESHOLD = 7
if MODIS:
ds = lib.data.modis.get_data(args, fixes=fixes).isel(time=0)
ds = ds.rename(chl_ocx='CHL')
else:
ds = lib.data.ostia.get_data(args, fixes=fixes).isel(time=0)
ds = ds.merge(lib.data.globcolour.get_data(
args, fixes=fixes)['CHL'].isel(time=0), compat='override')
ex = Example(-47.5, -45.8, 30.4, 32, (2007, 4, 7))
ds = ds.merge(lib.data.hi.get_data(args, fixes=fixes).isel(time=0))
ds = ex.extract(ds)
args = dict(region='GS', days=1,
scale=30., number=2, coef=0)
fixes = dict(Y=ex.date.year, m=ex.date.month, d=ex.date.day)
ds['HI'] = lib.data.hi.apply_coef(ds, lib.data.hi.get_coef(args))
sst = lo.get_data(args, fixes=fixes).sst.isel(time=0)
hi = lhi.get_data(args, fixes=fixes).isel(time=0)
im_kw = dict(add_labels=False, add_colorbar=False)
sst = ex.extract(sst)
hi = ex.extract(hi)
fig, axes = plt.subplots(2, 3, figsize=(5, 4.4),
sharex=True, sharey=True)
fig.subplots_adjust(left=0.08, bottom=0.05, right=0.98, top=0.90,
wspace=0.1, hspace=0.3)
coef = lhi.get_coef(args)
hi['HI'] = lhi.apply_coef(hi, coef)
hi['S'] = np.fabs(hi.S)
for v in 'SVB':
hi[v] = hi[v] * coef[v] * coef['HI']
ax1, ax2, ax3 = axes[0]
vmax = 22
c_lim = dict()
if ex_id == 0:
c_lim.update(dict(vmin=0.1, vmax=0.25))
im_kw = dict(add_labels=False, add_colorbar=False)
fig, axes = plt.subplots(2, 3, figsize=(5., 3.3), sharex=True, sharey=True)
fig.subplots_adjust(left=0.08, bottom=0.05, right=0.95, top=0.90,
wspace=0.3, hspace=0.3)
im1 = ds.sst.plot.imshow(ax=ax1, **im_kw, cmap=cm.thermal)
im2 = ds.CHL.plot.imshow(ax=ax2, **im_kw, **c_lim, cmap=cm.algae)
im3 = ds.HI.plot.imshow(ax=ax3, **im_kw, cmap=cm.haline)
im_sst = sst.plot.imshow(ax=axes[0, 0], **im_kw, cmap=cm.thermal)
images = [im1, im2, im3]
for i, var in enumerate(['S', 'V', 'B']):
images.append(ds[var].plot.imshow(ax=axes[1][i], **im_kw, cmap=cm.haline))
im_hi = hi.HI.plot.imshow(ax=axes[0, 1], **im_kw, cmap=cm.haline,
vmin=0, vmax=vmax)
labels = [r'SST [\si{\degreeCelsius}]',
r'Chl-\emph{a} [\(\si{\milli\gram\per\cubic\meter}\)]',
r'HI',
'Skewness', 'Variance', 'Bimodality']
for ax, im, label in zip(axes.flat, images, labels):
cax = plot_util.add_cax(ax, loc="top", pad=0.01, size=0.08)
fig.colorbar(im, cax=cax, ax=ax, orientation='horizontal', label=label,
ticklocation='top')
for ax, v in zip(axes[1], 'SVB'):
im_c = hi[v].plot.imshow(ax=ax, **im_kw, cmap=cm.haline,
vmin=0, vmax=vmax/3)
for ax in axes[0]:
ds['HI'].plot.contour(ax=ax, add_labels=False, add_colorbar=False,
levels=[THRESHOLD], colors='k')
labels = ['SST', 'HI', None, 'Skewness', 'Variance', 'Bimodality']
for ax, lab in zip(axes.flat, labels):
ax.set_title(lab)
for ax in axes.flat:
hi.HI.plot.contour(ax=ax, levels=[THRESHOLD], add_labels=False,
colors='k', linewidths=0.7)
ax.set_aspect('equal')
ax.xaxis.set_major_locator(MultipleLocator(1.))
ax.yaxis.set_major_locator(MultipleLocator(1.))
ax.xaxis.set_major_formatter(plot_util.LonFormatter('.0f'))
ax.yaxis.set_major_formatter(plot_util.LatFormatter('.0f'))
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])
cbar = fig.colorbar(im, cax=cax)
return cbar
# for ax in axes[0, :2]:
# win.add_to_plot(ax, ec='k', lw=0.7, fill=None, zorder=2.)
fig.canvas.draw()
add_cbar(im_sst, axes[0, 0])
add_cbar(im_hi, axes[0, 1])
add_cbar(im_c, axes[1, 2])
fig.savefig(path.join(lib.root_plot, 'Images', 'article_example_zoom',
'fig_ex{}_modis{:d}.png'.format(ex_id, MODIS)),
dpi=150)
'fig.png'), dpi=150)
Plots/article_climato.py
View file @ 5cac9485
......@@ -3,126 +3,56 @@ from os import path
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from cmocean import cm
from matplotlib.colors import LogNorm
from scipy import interpolate, ndimage
import lib
import lib.data.clouds_maps
import lib.data.ostia
import lib.data.globcolour
import lib.data.front_probability
import lib.data.p_frt_mask
import lib.zones
import lib.data.front_probability as lfrt
import lib.zones as lz
import Plots.util as plot_util
args = dict(region='GS')
sst = lib.data.ostia.get_data(args, climato=12).sst.mean('time')
chl = lib.data.globcolour.get_data(args, climato=12).CHL.mean('time')
frt = lib.data.front_probability.get_data(args, period='total').p_frt
clds = lib.data.clouds_maps.get_data(args) .cloud_probability.mean('time')
p_thr = lib.data.p_frt_mask.default_threshold
# land = lib.zones.get_land(args, grid='4km_EPSG32662').land
land = ~np.isfinite(frt) # Shady but works well eh
sst = sst.where(~land)
chl = chl.where(~land)
frt = frt.where(~land)
clds = clds.where(~land)
def filter(da, sigma):
def filter_ufunc(a, sigma=2):
x = np.arange(0, a.shape[1])
y = np.arange(0, a.shape[0])
xx, yy = np.meshgrid(x, y)
mask = np.isfinite(a)
x1, y1 = xx[mask], yy[mask]
fill = a[mask]
fill = interpolate.griddata((x1, y1), fill.ravel(), (xx, yy),
method='nearest')
filt = ndimage.gaussian_filter(fill, sigma)
filt[~mask] = np.nan
return filt
args = dict(region='GS', kind='2thr',
fixes=dict(thr_lo=5., thr_hi=10.))
da = xr.apply_ufunc(filter_ufunc, da,