Change argument to mask

Change histograms locations, add 'kind' directory
Fix all occurences

Compute/histogram/compute_hists_ostia_gs3.py

@@ -64,7 +64,7 @@ def main(args):
     m_next("Applying static masks")
     # Apply masks: land (enlarged), total zone, min front proba
     static = ~ds.land_large * ds.total
-    if not args['nomask']:
+    if args['mask']:
         static = static * ds.p_frt
     ds['HI'] = ds.HI.where(static)
     ds = ds.drop_vars(['land_large', 'total', 'p_frt'])
@@ -202,7 +202,7 @@ if __name__ == '__main__':
         parser.add_argument('-zones', type=str, default='INS')
 
     args = lib.get_args(['region', 'year', 'days', 'scale', 'number',
-                         'coef', 'fixes', 'threshold', 'nomask'], add_args)
+                         'coef', 'fixes', 'threshold', 'mask'], add_args)
     args['fixes']['Y'] = args['year']
     args['Y'] = args['year']
 

Compute/histogram/compute_hists_ostia_gs3_2thr.py

@@ -64,7 +64,7 @@ def main(args):
     m_next("Applying static masks")
     # Apply masks: land (enlarged), total zone, min front proba
     static = ~ds.land_large * ds.total
-    if not args['nomask']:
+    if args['mask']:
         static = static * ds.p_frt
     ds['HI'] = ds.HI.where(static)
     ds = ds.drop_vars(['land_large', 'total', 'p_frt'])
@@ -202,7 +202,7 @@ if __name__ == '__main__':
         parser.add_argument('-zones', type=str, default='INS')
 
     args = lib.get_args(['region', 'year', 'days', 'scale', 'number',
-                         'coef', 'fixes', 'thr_lo', 'thr_hi', 'nomask'],
+                         'coef', 'fixes', 'thr_lo', 'thr_hi', 'mask'],
                         add_args)
     args['fixes']['Y'] = args['year']
     args['Y'] = args['year']

Plots/Hists/all_zones_median.py

@@ -14,11 +14,9 @@ fixes = dict(
     scale=30.,
     number=2,
     coef=0,
-    thr_lo=6.,
+    thr_lo=5.,
     thr_hi=10.,
-    Y=2007,
-    nomask=False
-)
+    Y=2007)
 fixes.update(args['fixes'])
 args['fixes'] = fixes
 

Plots/Hists/all_zones_median_cli.py

@@ -10,14 +10,13 @@ import Plots.util as plot_util
 
 plot_util.use_tex(True)
 
-args = dict(region='GS', days=1)
+args = dict(region='GS', days=1, kind='2thr', mask=False)
 fixes = dict(
     scale=30.,
     number=2,
     coef=0,
-    thr_lo=6.,
-    thr_hi=10.,
-    nomask=False
+    thr_lo=5.,
+    thr_hi=10.
 )
 
 
@@ -57,10 +56,7 @@ def plot(cli, cli_d1, cli_d9, ax, var, sel, color='k', ls='-'):
                     color=color, alpha=.3)
 
 
-cset = tc.tol_cset('vibrant')
-colors = dict(low='red', mid='blue', hi='teal')
-
-colors = {m: getattr(cset, c) for m, c in colors.items()}
+colors = plot_util.get_mask_colors('vibrant', 'red', 'blue', 'teal')
 
 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,
@@ -98,10 +94,8 @@ axes[0, 0].set_title('South')
 axes[0, 1].set_title('Jet')
 axes[0, 2].set_title('North')
 
-labels = ['Low', 'Med', 'Hi']
-handles = [plt.Line2D((), (), ls='-', color=c) for c in colors.values()]
-fig.legend(handles, labels, bbox_to_anchor=[0.52, 0.78], loc='upper left',
-           framealpha=1.0)
+fig.legend(*plot_util.get_mask_legend(colors),
+           bbox_to_anchor=[0.52, 0.78], loc='upper left', framealpha=1.0)
 
 fig.canvas.draw()
 # fig.savefig(path.join(lib.root_plot, 'Hists', 'median_GS3', 'fig.pdf'))

Plots/Hists/article_median.py

 
-import matplotlib.pyplot as plt
-from matplotlib.dates import DateFormatter
-import pandas as pd
-import tol_colors as tc
+from os import path
 
+import matplotlib.pyplot as plt
 
 import lib
-import lib.data.hists
+import lib.data.hists as lh
+import Plots.util as plot_util
 
+plot_util.set_style()
+plot_util.use_tex(True)
 
-args = dict(
-    region='GS',
-    days=1
-)
+args = dict(region='GS', days=1, mask=True, kind='2thr')
 
+start = 2007
 fixes = dict(
     number=2,
     scale=30,
     coef=0,
-    zone='GS_S',
-    threshold=6.0
+    thr_lo=5.0,
+    thr_hi=10.,
+    Y=list(range(start, start+3))
 )
 
-variable = 'CHL'
-
-ds = lib.data.hists.get_data(args, fixes=fixes)
-ds = ds.sel(variable=variable)
+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')
+med = lh.get_median(pdf)
+med = med.load()
 
-ds.load()
-
-lib.data.hists.normalize_hist(ds)
-med = lib.data.hists.get_percentile(ds, 0.5)
-
-grp = med.groupby('time.dayofyear')
-cli = grp.mean()
-cli_std = grp.std()
-
-cli = cli.rename(dayofyear='time')
-cli_std = cli_std.rename(dayofyear='time')
+ymax = med.sel(mask='hi').quantile(0.99, 'time')
 
-cli = cli.assign_coords(time=pd.to_datetime(2000*1000 + cli.time.values,
-                                            format='%Y%j'))
-cli_std = cli_std.assign_coords(time=pd.to_datetime(2000*1000 + cli_std.time.values,
-                                                    format='%Y%j'))
+colors = plot_util.get_mask_colors('vibrant', 'red', 'blue', 'teal')
 
-# cli = cli.resample(time='8D').mean()
-# cli = cli.isel(time=slice(None, -1))
+fig, axes = plt.subplots(3, 1, figsize=(8, 4.5), sharex=True)
+right = 0.96
+fig.subplots_adjust(left=0.08, bottom=0.05, right=right, top=0.93,
+                    hspace=0.1)
 
-cset = tc.tol_cset('vibrant')
+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)
 
-fig, ax = plt.subplots(figsize=(5, 4))
-fig.subplots_adjust(left=0.15, bottom=0.08, right=0.98, top=0.98)
+    if z != 'N':
+        ax.set_ylim(None, ymax.sel(zone=zone))
+    ax.annotate(lh.get_GS3_title(z), (1., 0.5), xycoords='axes fraction',
+                xytext=(10, 0), textcoords='offset points',
+                ha='left', va='center', rotation='vertical')
 
+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')
+    ax.xaxis.axis_date()
+    ax.set_ylabel(plot_util.CHL_LABEL, labelpad=10.0)
 
-def plot(mask, color):
-    values = cli.sel(mask=mask)
-    v_std = cli_std.sel(mask=mask)
-    ax.plot(cli.time, values, color=color)
-    ax.fill_between(cli.time, values-2.*v_std, values+2.*v_std,
-                    color=color, alpha=0.4)
+axes[-1].set_xlim(*med.time[[0, -1]].values)
 
+fig.legend(*plot_util.get_mask_legend(colors),
+           bbox_to_anchor=[right, 0.99], loc='upper right', ncol=3,
+           borderaxespad=0.)
 
-plot('frt', cset.orange)
-plot('bkg', cset.teal)
 
-ax.set_xlim(*cli.time[[0, -1]].values)
-ax.xaxis.set_major_formatter(DateFormatter('%b'))
+fig.savefig(path.join(lib.root_plot, 'Hists', 'article_timeline',
+                      'fig.png'),
+            dpi=150)

Plots/Hists/article_surplus.py

@@ -6,7 +6,6 @@ import tol_colors as tc
 import xarray as xr
 
 import lib
-import lib.data.valid_pixels as lv
 import lib.data.hists as lh
 import Plots.util as plot_util
 
@@ -26,29 +25,47 @@ ds = ds.sum('zone')
 pdf = lh.get_pdf(ds, variables='CHL')
 med = lh.get_median(pdf).load()
 
-fixes['nomask'] = True
-ds = lh.get_data(fixes=fixes)
-ds = ds.squeeze()
-ds = lh.resample_8D_multiple_years(ds)
-ds = ds.sum('zone')
+chl = lh.get_data(fixes=fixes, mask=False).hist_CHL
+chl = chl.squeeze()
+chl = lh.resample_8D_multiple_years(chl)
+chl = chl.sum('zone', keep_attrs=True)
+
+T = lh.get_total(chl)
+N = lh.get_N(chl)
 
-T = lh.get_total(ds, variables='CHL')
-N = lh.get_N(ds, variables='CHL')
 
-T_f = T.sum('mask')
-T_nf_mid = (T.sel(mask=['low', 'hi']).sum('mask')
-            + med.sel(mask='low') * N.sel(mask='mid'))
-T_nf_hi = (T.sel(mask=['low', 'mid']).sum('mask')
-           + med.sel(mask='low') * N.sel(mask='hi'))
-T_nf = xr.concat([T_nf_mid, T_nf_hi], dim=xr.DataArray(['mid', 'hi'],
+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']))
-# T_nf = T.sel(mask='low') + med * N
+surplus = (C-C_nf)/C_nf * 100.
+
+# T_f = T.sum('mask')
+# T_nf_mid = (T.sel(mask=['low', 'hi']).sum('mask')
+#             + med.sel(mask='low') * N.sel(mask='mid'))
+# T_nf_hi = (T.sel(mask=['low', 'mid']).sum('mask')
+#            + med.sel(mask='low') * N.sel(mask='hi'))
+# T_nf = xr.concat([T_nf_mid, T_nf_hi], dim=xr.DataArray(['mid', 'hi'],
+#                                                        dims=['mask']))
+# # T_nf = T.sel(mask='low') + med * N
+
+# surplus = (T_f-T_nf)/T_f * 100.
+# surplus = surplus.sel(mask=['mid', 'hi'])
 
-surplus = (T_f-T_nf)/T_f * 100.
-surplus = surplus.load()
-# surplus = surplus.sel(mask=['mid', 'hi']).load()
+frac = N / N.sum('mask')
+est = frac * (med - med.sel(mask='low'))
+est = est.sel(mask=['mid', 'hi'])
 
-clis = lh.get_climatology(surplus)
+clis = lh.get_climatology(surplus.load())
+clis_est = lh.get_climatology(est.load())
 
 
 def plot(cli, cli_d1, cli_d9, ax, sel, color='k', ls='-'):
@@ -72,15 +89,16 @@ ax.xaxis.set_major_formatter(plot_util.MonthFormatter())
 
 for m, c in colors.items():
     plot(*clis, ax, dict(mask=m), color=c)
+    plot(*clis_est, ax, dict(mask=m), color=c)
 
 ax.set_ylabel('Surplus [%]')
 ax.set_xlim(clis[0].time[[0, -1]].values)
-ax.set_ylim(0, None)
+# ax.set_ylim(0, None)
 
 labels = list(colors.keys())
 handles = [plt.Line2D((), (), color=c) for c in colors.values()]
 fig.legend(handles, labels, bbox_to_anchor=[0., 1.], loc='upper left',
            ncol=2)
 
-fig.savefig(path.join(lib.root_plot, 'Hists', 'surplus', 'article_surplus.png'),
-            dpi=150)
+# fig.savefig(path.join(lib.root_plot, 'Hists', 'surplus', 'article_surplus.png'),
+#             dpi=150)

Plots/Hists/hists_seasons.py

@@ -14,9 +14,8 @@ fixes = dict(
     scale=30.,
     number=2,
     coef=0,
-    thr_lo=6.,
-    thr_hi=10.,
-    nomask=True
+    thr_lo=5.,
+    thr_hi=10.
 )
 fixes.update(args['fixes'])
 args['fixes'] = fixes

Plots/Hists/surplus_zones.py

@@ -6,7 +6,6 @@ import tol_colors as tc
 import xarray as xr
 
 import lib
-import lib.data.valid_pixels as lv
 import lib.data.hists as lh
 import Plots.util as plot_util
 
@@ -15,32 +14,41 @@ fixes = dict(
     number=2,
     scale=30.,
     coef=0,
-    thr_lo=6.,
+    thr_lo=5.,
     thr_hi=10.
 )
 
-ds = lh.get_data(fixes=fixes)
-ds = ds.squeeze()
-ds = lh.resample_8D_multiple_years(ds)
-pdf = lh.get_pdf(ds, variables='CHL')
+def sum_masks(da, masks):
+    return da.sel(mask=masks).sum('mask', keep_attrs=True)
+
+
+def concat(*args, masks=['mid', 'hi']):
+    assert len(args) == len(masks), 'ouch'
+    out = xr.concat(args, dim=xr.DataArray(masks, dims=['mask']))
+    return out
+
+
+hist = lh.get_data(fixes=fixes).hist_CHL.squeeze()
+# hist = hist.sum('zone', keep_attrs=True)
+hist = lh.resample_8D_multiple_years(hist)
+
+pdf_mid = sum_masks(hist, ['low', 'hi'])
+pdf_hi = sum_masks(hist, ['low', 'mid'])
+pdf = concat(pdf_mid, pdf_hi)
+
 med = lh.get_median(pdf).load()
 
-fixes['nomask'] = True
-ds = lh.get_data(fixes=fixes)
-ds = ds.squeeze()
-ds = lh.resample_8D_multiple_years(ds)
+hist = lh.get_data(fixes=fixes, mask=False).hist_CHL.squeeze()
+# hist = hist.sum('zone', keep_attrs=True)
+hist = lh.resample_8D_multiple_years(hist)
 
-T = lh.get_total(ds, variables='CHL')
-N = lh.get_N(ds, variables='CHL')
+T = lh.get_total(hist, variables='CHL')
+N = lh.get_N(hist, variables='CHL')
 
 T_f = T.sum('mask')
-T_nf_mid = (T.sel(mask=['low', 'hi']).sum('mask')
-            + med.sel(mask='low') * N.sel(mask='mid'))
-T_nf_hi = (T.sel(mask=['low', 'mid']).sum('mask')
-           + med.sel(mask='low') * N.sel(mask='hi'))
-T_nf = xr.concat([T_nf_mid, T_nf_hi], dim=xr.DataArray(['mid', 'hi'],
-                                                       dims=['mask']))
-# T_nf = T.sel(mask='low') + med * N
+T_nf_mid = sum_masks(T, ['low', 'hi']) + med.sel(mask='mid') * N.sel(mask='mid')
+N_nf_hi = sum_masks(T, ['low', 'mid']) + med.sel(mask='hi') * N.sel(mask='hi')
+T_nf = concat(T_nf_mid, N_nf_hi)
 
 surplus = (T_f-T_nf)/T_f * 100.
 surplus = surplus.load()
@@ -77,7 +85,6 @@ for ax, (z, t) in zip(axes, zones.items()):
         plot(*clis, ax, sel, color=c)
 
     ax.set_xlim(clis[0].time[[0, -1]].values)
-    ax.set_ylim(0, 27)
     ax.set_title(t)
 
 axes[0].set_ylabel('Surplus [%]')

lib/data/hists.py

@@ -11,28 +11,32 @@ import lib
 import lib.data
 
 
-ARGS = {'region', 'days', 'kind', 'nomask'}
-defaults = dict(kind='2thr', nomask=False)
+ARGS = {'region', 'days', 'kind', 'mask'}
+defaults = dict(kind='2thr', mask=True)
+
+"""
+kind:
+- '2thr': with 2 threshold, gives 3 histograms: low, mid and hi
+- '1thr': one threshold, gives 2 histograms: bkg and frt
+- '2d': 2d histogram (for Chl and HI)
+
+mask: If it should use the front presence probability mask, or not
+"""
 
 
 def PREGEX(args):
     pregex = ["number_%(number:fmt=d:rgx=%I)/",
               "scale_%(scale:fmt=.1f)/",
-              "coef_%(coef:fmt=d)/"]
+              "coef_%(coef:fmt=d)/"
+              "{}/".format(args['kind'])]
 
     if args['kind'] == '2thr':
-        pregex += ["hist_GS3",
-                   "_thr_%(thr_lo:fmt=.2f)_%(thr_hi:fmt=.2f)"]
-    elif args['kind'] == '2d':
-        pregex.append("hist_2D_GS3")
-    elif args['kind'] == 'thr':
-        pregex += ["hist_GS3",
-                   "_thr_%(threshold:fmt=.2f)"]
-    else:
-        raise ValueError("Unknown hist pregex type '{}'".format(args['dim']))
+        pregex.append("_thr_%(thr_lo:fmt=.2f)_%(thr_hi:fmt=.2f)")
+    elif args['kind'] == '1thr':
+        pregex.append("_thr_%(threshold:fmt=.2f)")
 
     pregex += ["_%(time:Y)",
-               "_nomask" if args['nomask'] else "",
+               "_nomask" if not args['mask'] else "",
                 ".nc"]
 
     return ''.join(pregex)
@@ -270,3 +274,12 @@ def get_climatology(ds):
     cli_d1 = restore_time(cli_d1)
     cli_d9 = restore_time(cli_d9)
     return cli, cli_d1, cli_d9
+
+
+def get_GS3_title(z: str):
+    d = dict(S='South', I='Jet', N='North')
+    return d[z]
+
+
+def get_GS3_name(z: str):
+    return 'GS3_{}'.format(z)