Small changes

unseen/eva.py

@@ -400,6 +400,7 @@ def gev_return_curve(
     n_bootstraps=1000,
     max_return_period=4,
     user_estimates=None,
+    max_shape_ratio=None,
 ):
     """Return x and y data for a GEV return period curve.
 
@@ -414,6 +415,9 @@ def gev_return_curve(
         The maximum return period is 10^{max_return_period}
     user_estimates: list, default None
         Initial estimates of the shape, loc and scale parameters
+    max_shape_ratio: float, optional
+        Maximum bootstrap shape parameter to full population shape parameter ratio (e.g. 6.0)
+        Useful for filtering bad fits to bootstrap samples
     """
 
     # GEV fit to data
@@ -438,7 +442,10 @@ def gev_return_curve(
         elif bootstrap_method == "non-parametric":
             boot_data = np.random.choice(data, size=data.shape, replace=True)
         boot_shape, boot_loc, boot_scale = fit_gev(boot_data, generate_estimates=True)
-
+        if max_shape_ratio:
+            shape_ratio = abs(boot_shape) / abs(shape)
+            if shape_ratio > max_shape_ratio:
+                continue
         boot_value = genextreme.isf(
             curve_probabilities, boot_shape, boot_loc, boot_scale
         )
@@ -486,6 +493,7 @@ def plot_gev_return_curve(
     ylim=None,
     text=False,
     user_estimates=None,
+    max_shape_ratio=None,
 ):
     """Plot a single return period curve.
 
@@ -509,6 +517,9 @@ def plot_gev_return_curve(
        Write the return period (and 95% CI) on the plot
     user_estimates: list, default None
         Initial estimates of the shape, loc and scale parameters
+    max_shape_ratio: float, optional
+        Maximum bootstrap shape parameter to full population shape parameter ratio (e.g. 6.0)
+        Useful for filtering bad fits to bootstrap samples
     """
 
     if direction == "deceedance":
@@ -521,6 +532,7 @@ def plot_gev_return_curve(
         n_bootstraps=n_bootstraps,
         max_return_period=max_return_period,
         user_estimates=user_estimates,
+        max_shape_ratio=max_shape_ratio,
     )
     (
         curve_return_periods,

unseen/moments.py

@@ -3,6 +3,7 @@
 import argparse
 import logging
 
+import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 import scipy
@@ -12,6 +13,11 @@
 
 
 logging.basicConfig(level=logging.INFO)
+mpl.rcParams["axes.labelsize"] = "x-large"
+mpl.rcParams["axes.titlesize"] = "xx-large"
+mpl.rcParams["xtick.labelsize"] = "x-large"
+mpl.rcParams["ytick.labelsize"] = "x-large"
+mpl.rcParams["legend.fontsize"] = "x-large"
 
 
 def calc_ci(data):
@@ -178,7 +184,7 @@ def create_plot(
         "GEV scale": "scale parameter",
         "GEV location": "location parameter",
     }
-    fig = plt.figure(figsize=[15, 22])
+    fig = plt.figure(figsize=[15, 28])
     for plotnum, moment in enumerate(moments):
         ax = fig.add_subplot(4, 2, plotnum + 1)
         ax.hist(
@@ -217,8 +223,8 @@ def create_plot(
                 linestyle="--",
                 linewidth=3.0,
             )
-        ax.set_ylabel("count")
-        ax.set_xlabel(units[moment])
+        ax.set_ylabel("count", fontsize="large")
+        ax.set_xlabel(units[moment], fontsize="large")
         letter = letters[plotnum]
         ax.set_title(f"({letter}) {moment}")
         if letter == "a":
@@ -325,9 +331,9 @@ def _main():
         infile_logs = None
 
     create_plot(
-        da_fcst,
-        da_obs,
-        da_bc_fcst=da_bc_fcst,
+        da_fcst.compute(),
+        da_obs.compute(),
+        da_bc_fcst=da_bc_fcst.compute(),
         outfile=args.outfile,
         units=args.units,
         ensemble_dim=args.ensemble_dim,

unseen/stability.py

@@ -31,7 +31,7 @@ def plot_dist_by_lead(ax, sample_da, metric, units=None, lead_dim="lead_time"):
     """
 
     lead_times = np.unique(sample_da[lead_dim].values)
-    colors = iter(plt.cm.viridis_r(np.linspace(0, 1, len(lead_times))))
+    colors = iter(plt.cm.BuPu(np.linspace(0, 1, len(lead_times))))
     for lead in lead_times:
         selection_da = sample_da.sel({lead_dim: lead})
         selection_da = selection_da.dropna("sample")
@@ -156,7 +156,7 @@ def plot_return_by_lead(
     """
 
     lead_times = np.unique(sample_da["lead_time"].values)
-    colors = iter(plt.cm.viridis_r(np.linspace(0, 1, len(lead_times))))
+    colors = iter(plt.cm.BuPu(np.linspace(0, 1, len(lead_times))))
     for lead in lead_times:
         selection_da = sample_da.sel({"lead_time": lead})
         selection_da = selection_da.dropna("sample")
@@ -188,9 +188,9 @@ def plot_return_by_lead(
             alpha=0.3,
         )
 
-    ax.grid(True)
     ax.set_title(f"(b) {metric} return period by lead time")
     ax.set_xscale("log")
+    ax.grid(True, which="both")
     ax.set_xlabel("return period (years)")
     units_label = units if units else sample_da.attrs["units"]
     ax.set_ylabel(units_label)
@@ -260,9 +260,9 @@ def plot_return_by_time(
             alpha=0.3,
         )
 
-    ax.grid(True)
     ax.set_title(f"(d) {metric} return period by year")
     ax.set_xscale("log")
+    ax.grid(True, which="both")
     ax.set_xlabel("return period (years)")
     units_label = units if units else sample_da.attrs["units"]
     ax.set_ylabel(units_label)
@@ -428,7 +428,7 @@ def _main():
     da_fcst = ds_fcst[args.var]
 
     create_plot(
-        da_fcst,
+        da_fcst.compute(),
         args.metric,
         args.start_years,
         outfile=args.outfile,