Save resampled climates from KCS diagnostic local_resampling.py (#2221)

* store resampled climates to nc files

* minor style fixes

* added step in main to save resampled nc files

* added docstring

* fixed line length?

* add myself

* spelling mistake corrected

Co-authored-by: Manuel Schlund <manuel.schlund@dlr.de>

doc/sphinx/source/recipes/recipe_kcs.rst

@@ -148,6 +148,9 @@ The diagnostic ``local_resampling`` procudes a number of output files:
      Combination 7          6         3         2         1         6         1         2         1         0         2         1         3
 
 
+* ``resampled_control_<scenario>.nc``: containing the monthly means for the control period according to the final combinations.
+* ``resampled_future_<scenario>.nc``: containing the monthly means for the future period according to the final combinations.
+
 * Provenance information: bibtex, xml, and/or text files containing citation information are stored alongside the final result and the final figure.
   The final combinations only derive from the target model data, whereas the figure also uses CMIP data.
 * A figure used to validate the final result, reproducing figures 5 and 6 from Lenderink et al.:

esmvaltool/config-references.yml

@@ -154,6 +154,10 @@ authors:
     name: Mohr, Christian Wilhelm
     institute: Cicero, Norway
     orcid: https://orcid.org/0000-0003-2656-1802
+  daniels_emma:
+    name: Daniels, Emma
+    institute: KNMI, Netherlands
+    orcid:
   davin_edouardleopold:
     name: Davin, Edouard Leopold
     institute: ETH Zurich, Switzerland

esmvaltool/diag_scripts/kcs/local_resampling.py

@@ -28,6 +28,7 @@ def _create_provenance_record(ancestor_files):
             'kalverla_peter',
             'alidoost_sarah',
             'rol_evert',
+            'daniels_emma',
         ],
         'ancestors': ancestor_files,
     }
@@ -366,13 +367,17 @@ def _recombine(segments, combinations):
                                coords={'segment': range(n_segments)})
 
         # Recombine the segments using the indexer
-        resample = segments.sel(ensemble_member=indexer).mean('segment')
+        resample = segments.sel(ensemble_member=indexer).mean('segment',
+                                                              keep_attrs=True)
         new_climates.append(resample)
     return xr.concat(new_climates, dim='sample')
 
 
-def _get_climatology(cfg, scenario_name, table):
-    """Determine the change in <variable> PDF of each scenario."""
+def _get_climatology(cfg, scenario_name, table, prov=None):
+    """Determine the change in <variable> PDF of each scenario.
+
+    Save the resampled climates of each scenario to nc files.
+    """
     dataset, _ = _get_data_target_model(cfg)
 
     future = cfg['scenarios'][scenario_name]['resampling_period']
@@ -381,6 +386,29 @@ def _get_climatology(cfg, scenario_name, table):
 
     resampled_control = _recombine(segments_control, table['control'])
     resampled_future = _recombine(segments_future, table['future'])
+    # Store the resampled contol climates
+    filename = get_diagnostic_filename(f'resampled_control_{scenario_name}',
+                                       cfg,
+                                       extension='nc')
+    resampled_control.to_netcdf(filename)
+    LOGGER.info("Created control resamples for scenario %s: \n %s",
+                scenario_name, table)
+    LOGGER.info('Output stored as %s', filename)
+    # # Write provenance information
+    with ProvenanceLogger(cfg) as provenance_logger:
+        provenance_logger.log(filename, prov)
+
+        # Store the resampled future climates
+    filename = get_diagnostic_filename(f'resampled_future_{scenario_name}',
+                                       cfg,
+                                       extension='nc')
+    resampled_future.to_netcdf(filename)
+    LOGGER.info("Created future resamples for scenario %s: \n %s",
+                scenario_name, table)
+    LOGGER.info('Output stored as %s', filename)
+    # # Write provenance information
+    with ProvenanceLogger(cfg) as provenance_logger:
+        provenance_logger.log(filename, prov)
 
     quantiles = [.05, .1, .25, .5, .75, .90, .95]
     qcontrol = resampled_control.groupby('time.season').quantile(
@@ -393,17 +421,22 @@ def _get_climatology(cfg, scenario_name, table):
     return xr.merge([qchange_tas, qchange_pr])
 
 
-def make_plots(cfg, scenario_tables):
+def get_climatologies(cfg, scenario_tables, prov=None):
+    """Determine the changes in <variable> PDF of all scenarios."""
+    climates = {}
+    for name in cfg['scenarios'].keys():
+        climatology = _get_climatology(cfg, name, table=scenario_tables[name],
+                                       prov=prov)
+        climates[name] = climatology
+    return climates
+
+
+def make_plots(cfg, climates):
     """Reproduce figure 5 from the paper."""
     # Note that quantile is applied twice! Once to get the pdf's of seasonal
     # tas/pr and once to get the multimodel pdf of the quantile changes
     metadata = cfg['input_data'].values()
 
-    climates = {}
-    for name, info in cfg['scenarios'].items():
-        climatology = _get_climatology(cfg, name, table=scenario_tables[name])
-        climates[name] = climatology
-
     for year in [2050, 2085]:
         fig, subplots = plt.subplots(2, 2, figsize=(12, 8))
 
@@ -453,10 +486,13 @@ def main(cfg):
     subsets = get_percentile_subsets(cfg, segment_season_means, top1000s)
 
     # Step 3: select final set of eight samples
-    scenarios = select_final_subset(cfg, subsets, prov=provenance)
+    scenario_tables = select_final_subset(cfg, subsets, prov=provenance)
+
+    # Step 4: create the resampled climates
+    climates = get_climatologies(cfg, scenario_tables, prov=provenance)
 
-    # Step 4: plot the results
-    make_plots(cfg, scenarios)
+    # Step 5: plot the results
+    make_plots(cfg, climates)
 
 
 if __name__ == '__main__':