Merge pull request #66 from DamienIrving/master

Big merge of my fork

file_lists/CAFE_c5-d60-pX-f6_config.mk

@@ -0,0 +1,12 @@
+# Configuration file for Makefile workflows
+
+MODEL=CAFE
+EXPERIMENT=c5-d60-pX-f6
+BASE_PERIOD=1995-01-01 2020-12-31
+BASE_PERIOD_TEXT=1995-2020
+TIME_PERIOD_TEXT=19950501-20201101
+STABILITY_START_YEARS=1995 2000 2005 2010 2015
+MODEL_IO_OPTIONS=--metadata_file /home/599/dbi599/unseen/config/dataset_cafe_daily.yml
+MODEL_NINO_OPTIONS=--lon_bnds -170 -120 --lat_dim geolat_t --lon_dim geolon_t --agg_x_dim xt_ocean --agg_y_dim yt_ocean --anomaly ${BASE_PERIOD} --anomaly_freq month --metadata_file /home/599/dbi599/unseen/config/dataset_cafe_monthly_ocean.yml
+
+

file_lists/CMCC-CM2-SR5_dcppA-hindcast_config.mk

@@ -0,0 +1,12 @@
+# Configuration file for Makefile workflows
+
+MODEL=CMCC-CM2-SR5
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1960-01-01 2019-12-31
+BASE_PERIOD_TEXT=196011-201911
+TIME_PERIOD_TEXT=196011-201911
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 10
+MODEL_NINO_OPTIONS=${MODEL_IO_OPTIONS} --lon_bnds 190 240 --lat_dim latitude --lon_dim longitude --agg_y_dim j --agg_x_dim i --anomaly ${BASE_PERIOD} --anomaly_freq month
+
+

file_lists/CanESM5_dcppA-hindcast_config.mk

@@ -0,0 +1,14 @@
+# Configuration file for Makefile workflows
+
+MODEL=CanESM5
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1961-01-01 2017-12-31
+BASE_PERIOD_TEXT=196101-201701
+TIME_PERIOD_TEXT=196101-201701
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 20
+MODEL_NINO_OPTIONS=${MODEL_IO_OPTIONS} --lon_bnds 190 240 --lat_dim latitude --lon_dim longitude --agg_y_dim j --agg_x_dim i --anomaly ${BASE_PERIOD} --anomaly_freq month
+
+
+
+

file_lists/EC-Earth3_dcppA-hindcast_config.mk

@@ -0,0 +1,12 @@
+# Configuration file for Makefile workflows
+
+MODEL=EC-Earth3
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1960-01-01 2017-12-20
+BASE_PERIOD_TEXT=1960-2017
+TIME_PERIOD_TEXT=196011-201711
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 15 --n_time_files 11 --standard_calendar
+MODEL_NINO_OPTIONS=${MODEL_IO_OPTIONS} --lon_bnds 190 240 --lat_dim latitude --lon_dim longitude --agg_y_dim j --agg_x_dim i --anomaly ${BASE_PERIOD} --anomaly_freq month
+
+

file_lists/HadGEM3-GC31-MM_dcppA-hindcast_config.mk

@@ -0,0 +1,12 @@
+# Configuration file for Makefile workflows
+
+MODEL=HadGEM3-GC31-MM
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1970-01-01 2018-12-20
+BASE_PERIOD_TEXT=197011-201811
+TIME_PERIOD_TEXT=196011-201811
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 10 --n_time_files 12
+MODEL_NINO_OPTIONS=${MODEL_IO_OPTIONS} --lon_bnds -170 -120 --lat_dim latitude --lon_dim longitude --agg_y_dim j --agg_x_dim i --anomaly ${BASE_PERIOD} --anomaly_freq month
+
+

file_lists/IPSL-CM6A-LR_dcppA-hindcast_config.mk

@@ -0,0 +1,11 @@
+# Configuration file for Makefile workflows
+
+MODEL=IPSL-CM6A-LR
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1970-01-01 2017-12-31
+BASE_PERIOD_TEXT=197001-201701
+TIME_PERIOD_TEXT=196101-201701
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 10
+MODEL_NINO_OPTIONS=--n_ensemble_files 10 --lon_bnds -170 -120 --lat_dim nav_lat --lon_dim nav_lon --agg_y_dim y --agg_x_dim x --anomaly ${BASE_PERIOD} --anomaly_freq month
+

file_lists/MIROC6_dcppA-hindcast_config.mk

@@ -0,0 +1,14 @@
+
+MODEL=MIROC6
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1960-01-01 2018-12-31
+BASE_PERIOD_TEXT=196011-201811
+TIME_PERIOD_TEXT=196011-201811
+MIN_LEAD=1
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 10 --n_time_files 2
+MODEL_NINO_OPTIONS=--n_ensemble_files 10 --lon_bnds 190 240 --lat_dim latitude --lon_dim longitude --agg_y_dim y --agg_x_dim x --anomaly ${BASE_PERIOD} --anomaly_freq month
+
+
+
+

file_lists/MPI-ESM1-2-HR_dcppA-hindcast_config.mk

@@ -0,0 +1,11 @@
+# Configuration file for Makefile workflows
+
+MODEL=MPI-ESM1-2-HR
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1960-01-01 2018-12-31
+BASE_PERIOD_TEXT=196011-201811
+TIME_PERIOD_TEXT=196011-201811
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 10
+
+

file_lists/MRI-ESM2-0_dcppA-hindcast_config.mk

@@ -0,0 +1,14 @@
+# Configuration files for Makefile workflows
+
+MODEL=MRI-ESM2-0
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1965-01-01 2018-12-31
+BASE_PERIOD_TEXT=196011-201911
+TIME_PERIOD_TEXT=196011-201911
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 10
+MODEL_NINO_OPTIONS=--n_ensemble_files 10 --lon_bnds 190 240 --lat_dim latitude --lon_dim longitude --agg_y_dim y --agg_x_dim x --anomaly ${BASE_PERIOD} --anomaly_freq month
+
+
+
+

file_lists/NorCPM1_dcppA-hindcast_config.mk

@@ -0,0 +1,12 @@
+# Configuration file for Makefile workflows
+
+MODEL=NorCPM1
+EXPERIMENT=dcppA-hindcast
+BASE_PERIOD=1960-01-01 2018-12-31
+BASE_PERIOD_TEXT=1960-2018
+TIME_PERIOD_TEXT=196010-201810
+STABILITY_START_YEARS=1960 1970 1980 1990 2000 2010
+MODEL_IO_OPTIONS=--n_ensemble_files 20
+MODEL_NINO_OPTIONS=${MODEL_IO_OPTIONS} --lon_bnds 190 240 --lat_dim latitude --lon_dim longitude --agg_y_dim j --agg_x_dim i --anomaly ${BASE_PERIOD} --anomaly_freq month
+
+

unseen/bias_correction.py

@@ -227,7 +227,11 @@ def _main():
 
     if args.output_chunks:
         ds_fcst_bc = ds_fcst_bc.chunk(args.output_chunks)
-    fileio.to_zarr(ds_fcst_bc, args.outfile)
+
+    if "zarr" in args.outfile:
+        fileio.to_zarr(ds_fcst_bc, args.outfile)
+    else:
+        ds_fcst_bc.to_netcdf(args.outfile)
 
 
 if __name__ == "__main__":

unseen/eva.py

@@ -601,6 +601,8 @@ def gev_return_curve(
     bootstrap_method="non-parametric",
     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.
 
@@ -609,14 +611,22 @@ def gev_return_curve(
     data : xarray DataArray
     event_value : float
         Magnitude of event of interest
-    bootstrap_method : {'parametric', 'non-parametric'}, default 'non-parametric'
+    bootstrap_method : {'parametric', 'non-parametric'}, default "non-parametric"
     n_bootstraps : int, default 1000
     max_return_period : float, default 4
         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
-    shape, loc, scale = fit_gev(data, generate_estimates=True, stationary=True)
+    if user_estimates:
+        shape, loc, scale = fit_gev(data, user_estimates=user_estimates, stationary=True)
+    else:
+        shape, loc, scale = fit_gev(data, generate_estimates=True, stationary=True)
 
     curve_return_periods = np.logspace(0, max_return_period, num=10000)
     curve_probabilities = 1.0 / curve_return_periods
@@ -634,7 +644,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
         )
@@ -681,6 +694,8 @@ def plot_gev_return_curve(
     ylabel=None,
     ylim=None,
     text=False,
+    user_estimates=None,
+    max_shape_ratio=None,
 ):
     """Plot a single return period curve.
 
@@ -702,6 +717,11 @@ def plot_gev_return_curve(
         Limits for y-axis
     text : bool, default False
        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":
@@ -713,6 +733,8 @@ def plot_gev_return_curve(
         bootstrap_method=bootstrap_method,
         n_bootstraps=n_bootstraps,
         max_return_period=max_return_period,
+        user_estimates=user_estimates,
+        max_shape_ratio=max_shape_ratio,
     )
     (
         curve_return_periods,

unseen/fileio.py

@@ -46,10 +46,10 @@ def open_dataset(
     time_freq=None,
     time_agg=None,
     time_agg_dates=False,
-    month=None,
+    months=None,
     season=None,
     reset_times=False,
-    complete_time_agg_periods=False,
+    time_agg_min_tsteps=None,
     input_freq=None,
     time_dim="time",
     isel={},
@@ -111,14 +111,14 @@ def open_dataset(
         Record the date of each time aggregated event (e.g. annual max)
     standard_calendar : bool, default False
         Force a common calendar on all input files
-    month : {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}, optional
-        Select a single month from the dataset
+    months : list, optional
+        Select months from the dataset
     season : {'DJF', 'MAM', 'JJA', 'SON'}, optional
         Select a single season after Q-NOV temporal resampling
     reset_times : bool, default False
         Shift time values after resampling so months match initial date
-    complete_time_agg_periods : bool default False
-        Limit temporal aggregation output to complete years/months
+    time_agg_min_tsteps : int, optional
+        Minimum number of timesteps for temporal aggregation
     input_freq : {'A', 'Q', 'M', 'D'}, optional
         Input time frequency for resampling (estimated if not provided)
     time_dim: str, default 'time'
@@ -161,9 +161,9 @@ def open_dataset(
         ds = ds.isel(isel)
     if sel:
         ds = ds.sel(sel)
-    if month:
-        ds = time_utils.select_month(
-            ds, month, init_month=reset_times, time_dim=time_dim
+    if months:
+        ds = time_utils.select_months(
+            ds, months, init_month=reset_times, time_dim=time_dim
         )
 
     # Scale factors
@@ -231,7 +231,7 @@ def open_dataset(
             variables,
             season=season,
             reset_times=reset_times,
-            complete=complete_time_agg_periods,
+            min_tsteps=time_agg_min_tsteps,
             agg_dates=time_agg_dates,
         )
 
@@ -628,11 +628,11 @@ def _parse_command_line():
         help="Record the date of each time aggregated event (e.g. annual max) [default=False]",
     )
     parser.add_argument(
-        "--month",
+        "--months",
         type=int,
-        choices=(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12),
+        nargs='*',
         default=None,
-        help="Select a single month from the dataset",
+        help="Select months from the dataset",
     )
     parser.add_argument(
         "--season",
@@ -648,10 +648,10 @@ def _parse_command_line():
         help="Shift time values after resampling so months match initial date [default=False]",
     )
     parser.add_argument(
-        "--complete_time_agg_periods",
-        action="store_true",
-        default=False,
-        help="Limit temporal aggregation output to complete years/months [default=False]",
+        "--time_agg_min_tsteps",
+        type=int,
+        default=None,
+        help="Minimum number of timesteps for temporal aggregation [default=None]",
     )
     parser.add_argument(
         "--input_freq",
@@ -844,10 +844,10 @@ def _main():
         "time_freq": args.time_freq,
         "time_agg": args.time_agg,
         "time_agg_dates": args.time_agg_dates,
-        "month": args.month,
+        "months": args.months,
         "season": args.season,
         "reset_times": args.reset_times,
-        "complete_time_agg_periods": args.complete_time_agg_periods,
+        "time_agg_min_tsteps": args.time_agg_min_tsteps,
         "input_freq": args.input_freq,
         "time_dim": args.time_dim,
         "isel": args.isel,

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/similarity.py

@@ -216,8 +216,8 @@ def _main():
     ds_fcst = fileio.open_dataset(args.fcst_file, variables=[args.var])
     ds_obs = fileio.open_dataset(args.obs_file, variables=[args.var])
     if args.reference_time_period:
-        time_slice = time_utils.date_pair_to_time_slice(args.reference_time_period)
-        ds_obs = ds_obs.sel({args.time_dim: time_slice})
+        start_date, end_date = args.reference_time_period
+        ds_obs = ds_obs.sel({args.time_dim: slice(start_date, end_date)})
 
     ds_similarity = similarity_tests(
         ds_fcst,
@@ -239,7 +239,11 @@ def _main():
 
     if args.output_chunks:
         ds_similarity = ds_similarity.chunk(args.output_chunks)
-    fileio.to_zarr(ds_similarity, args.outfile)
+
+    if "zarr" in args.outfile:
+        fileio.to_zarr(ds_similarity, args.outfile)
+    else:
+        ds_similarity.to_netcdf(args.outfile)
 
 
 if __name__ == "__main__":