Diagnostic for calculating Lamb Weathertypes

doc/sphinx/source/recipes/index.rst

@@ -64,6 +64,7 @@ Atmosphere
    recipe_validation
    recipe_radiation_budget
    recipe_aod_aeronet_assess
+   recipe_weathertyping
 
 Climate metrics
 ^^^^^^^^^^^^^^^

doc/sphinx/source/recipes/recipe_weathertyping.rst

@@ -0,0 +1,112 @@
+.. _recipes_weathertyping:
+
+Lamb Weathertypes
+=====
+
+Overview
+--------
+
+A diagnostic to calculate Lamb weathertypes over a given region. Furthermore, 
+correlations between weathertypes and precipitation patterns over a given are can be calculated 
+and 'combined' or 'simplified' weathertypes can be derived. Additionally, mean fields, as well as
+anomalies and standard deviations can be plotted.
+
+
+Available recipes and diagnostics
+---------------------------------
+
+Recipes are stored in esmvaltool/recipes/
+
+    * recipe_weathertyping.yml
+
+Diagnostics are stored in esmvaltool/diag_scripts/weathertyping/
+
+    * weathertyping.py: calculate lamb and simplified WT, plot mean, anomalies and std for each WT for psl, tas, prcp
+
+
+User settings in recipe
+-----------------------
+
+#. weathertyping.py
+
+   *Required settings for script*
+
+   *Optional settings for script*
+
+   * correlation_threshold: correlation threshold for selecting similar WT pairs, only needed if automatic_slwt==True and predefined_slwt==False. default=0.9
+   * rmse_threshold: rmse threshold for selecting similar WT pairs, only needed if automatic_slwt==True and predefined_slwt==False. default=0.002
+   * plotting: if true, create plots of means, anomalies and std for psl, tas, prcp
+   * automatic_slwt: if true, automatically combine WT with similar precipitation patterns over specified area (via thresholds of correlation and rmse OR via predefined_slwt)
+   * predefined_slwt: dictionary of mappings between weathertypes
+
+.. note::
+
+  predefined_slwt can be a dictionary where keys are slwt and the values are arrays of lwt OR where keys are lwt and values are slwt
+
+   *Required settings for variables*
+
+   *Optional settings for variables*
+
+   *Required settings for preprocessor*
+
+   *Optional settings for preprocessor*
+
+   *Color tables*
+
+
+Variables
+---------
+
+* psl (atmos, day, time longitude latitude)
+* tas (atmos, day, time longitude latitude)
+* tp (atmos, day, time longitude latitude)
+* pr (atmos, day, time longitude latitude)
+
+
+Observations and reformat scripts
+---------------------------------
+
+*Note: (1) obs4MIPs data can be used directly without any preprocessing;
+(2) see headers of reformat scripts for non-obs4MIPs data for download
+instructions.*
+
+* E-OBS: European Climate Assessment & Dataset gridded observationl daily precipitation sum
+
+References
+----------
+
+* Maraun, D., Truhetz, H., & Schaffer, A. (2021). Regional climate model biases, their dependence on synoptic circulation biases and the potential for bias adjustment: A process-oriented evaluation of the Austrian regional climate projections. Journal of Geophysical Research: Atmospheres, 126, e2020JD032824. https://doi.org/10.1029/2020JD032824 
+* Jones, P.D., Hulme, M. and Briffa, K.R. (1993), A comparison of Lamb circulation types with an objective classification scheme. Int. J. Climatol., 13: 655-663. https://doi.org/10.1002/joc.3370130606
+
+Example plots
+-------------
+
+.. _fig_weathertyping_1:
+.. figure::  /recipes/figures/weathertyping/lwt_1_ERA5_psl_mean.png
+   :align:   center
+
+   PSL mean map of Lamb WT 1 for ERA5.
+
+.. _fig_weathertyping_2:
+.. figure::  /recipes/figures/weathertyping/lwt_1_TaiESM1_psl_mean.png
+   :align:   center
+
+   PSL mean map of Lamb WT 1 for TaiESM1.
+
+.. _fig_weathertyping_3:
+.. figure::  /recipes/figures/weathertyping/slwt_EOBS_4_ERA5_psl_mean.png
+   :align:   center
+
+   PSL mean map of slwt_EOBS 4 for ERA5 (in this case combined Lamb WT 24 and 23).
+
+.. _fig_weathertyping_4:
+.. figure::  /recipes/figures/weathertyping/correlation_matrix_E-OBS.png
+   :align:   center
+
+   Heatmap of correlation values for Lamb WTs 1-27.
+
+.. _fig_weathertyping_5:
+.. figure::  /recipes/figures/weathertyping/ERA5__lwt_rel_occurrence.png
+   :align:   center
+
+   Stackplot of seasonal relative occureences of each WT for ERA5.

esmvaltool/config-references.yml

@@ -286,6 +286,11 @@ authors:
     name: Juckes, Martin
     institute: BADC, UK
     orcid:
+  jury_martin:
+    name: Jury, Martin
+    institute: WEGC, Austria
+    e-mail: [email protected]
+    orcid: https://orcid.org/0000-0003-0590-7843
   kadygrov_nikolay:
     name: Kadygrov, Nikolay
     institute: IPSL, France
@@ -315,6 +320,12 @@ authors:
     name: Krasting, John
     institute: NOAA, USA
     orcid: https://orcid.org/0000-0002-4650-9844
+  kroissenbrunner_thomas:
+    name: Kroissenbrunner, Thomas
+    institute: WEGC, Austria
+    email: [email protected]
+    github: thomaskroi1996
+    orcid: 
   kuehbacher_birgit:
     name: Kuehbacher, Birgit
     institute: DLR, Germany
@@ -793,6 +804,7 @@ projects:
   trr181: DFG Project TRR-181, Energy transfers in Atmosphere and Ocean
   ukesm: UKESM, UK Earth System Model project (NERC)
   usmile: ERC Synergy Grant USMILE
+  preval: PREVAL ÖKS
 
 
 realms:

esmvaltool/diag_scripts/weathertyping/weathertyping.py

@@ -0,0 +1,222 @@
+from wt_utils import (get_cfg_vars, load_wt_preprocessors,
+                      wt_algorithm, get_ancestors_era5_eobs, calc_slwt_obs,
+                      run_predefined_slwt, combine_wt_to_file, load_wt_files,
+                      get_looping_dict, plot_means, get_model_output_filepath,
+                      calc_lwt_slwt_model, write_lwt_to_file, calc_lwt_model,
+                      plot_seasonal_occurrence
+                      )
+
+import iris
+
+# import internal esmvaltool modules here
+from esmvaltool.diag_scripts.shared import (
+    run_diagnostic,
+)
+
+
+def run_automatic_slwt(cfg: dict):
+    preproc_variables_dict, _, _, \
+        work_dir, plotting, _, predefined_slwt = get_cfg_vars(cfg)
+    for dataset_name, dataset_vars in preproc_variables_dict.items():
+        if dataset_name == 'ERA5':
+            wt_preproc, wt_preproc_prcp, wt_preproc_prcp_eobs = \
+                load_wt_preprocessors(dataset_name, preproc_variables_dict)
+
+            # calculate lwt
+            lwt = wt_algorithm(wt_preproc, dataset_name)
+
+            era5_ancestors, eobs_ancestors = get_ancestors_era5_eobs(
+                dataset_name, preproc_variables_dict)
+
+            # calculate simplified lwt based on precipitation
+            # patterns or use predefined_slwt
+            if predefined_slwt is False:
+                slwt_era5 = calc_slwt_obs(
+                    cfg,
+                    lwt,
+                    wt_preproc_prcp,
+                    dataset_name,
+                    era5_ancestors,
+                )
+                slwt_eobs = calc_slwt_obs(
+                    cfg,
+                    lwt,
+                    wt_preproc_prcp_eobs,
+                    'E-OBS',
+                    eobs_ancestors,
+                )
+            else:
+                slwt_era5, slwt_eobs = run_predefined_slwt(work_dir,
+                                                           dataset_name,
+                                                           lwt,
+                                                           predefined_slwt)
+
+            # write to file
+            wt_list = [lwt, slwt_era5, slwt_eobs]
+            combine_wt_to_file(cfg, wt_list, wt_preproc,
+                               dataset_name)
+
+            # load weathertype files as cubes
+            wt_cubes = load_wt_files(f'{work_dir}/{dataset_name}.nc')
+
+            var_dict = get_looping_dict(
+                dataset_vars
+            )  # dataset_vars is list of variables for dataset dataset_name
+            if plotting:
+                # plot means
+                for var_name, var_data in var_dict.items():
+                    data_info = {'dataset': dataset_name,
+                                 'var': var_name,
+                                 'preproc_path': var_data[1]}
+                    plot_means(cfg, var_data[0], wt_cubes, data_info)
+                plot_seasonal_occurrence(cfg, wt_cubes, dataset_name)
+        else:
+            if dataset_name == 'E-OBS':
+                continue
+            wt_preproc = iris.load_cube(dataset_vars[0].get('filename'))
+
+            output_file_path, preproc_path = get_model_output_filepath(
+                dataset_name, dataset_vars)
+
+            # calculate weathertypes
+            data_info = {'dataset': dataset_name,
+                         'preproc_path': preproc_path,
+                         'output_file_path': output_file_path}
+            calc_lwt_slwt_model(cfg, wt_preproc, data_info, predefined_slwt)
+
+            # load wt files
+            wt_cubes = load_wt_files(
+                f'{work_dir}/{output_file_path}/{dataset_name}.nc')
+
+            var_dict = get_looping_dict(
+                dataset_vars
+            )  # dataset_vars is list of variables for dataset dataset_name
+
+            if plotting:
+                # plot means
+                ensemble = dataset_vars[0].get('ensemble')
+                for var_name, var_data in var_dict.items():
+                    data_info = {'dataset': dataset_name,
+                                 'var': var_name,
+                                 'preproc_path': var_data[1],
+                                 'ensemble': ensemble}
+                    plot_means(cfg, var_data[0], wt_cubes, data_info,
+                               only_lwt=True)
+                plot_seasonal_occurrence(cfg, wt_cubes, dataset_name,
+                                         ensemble=ensemble)
+
+
+def run_lwt(cfg: dict):
+    preproc_variables_dict, _, _, \
+        work_dir, plotting, _, _ = get_cfg_vars(cfg)
+    for dataset_name, dataset_vars in preproc_variables_dict.items():
+        if dataset_name == 'ERA5':
+            wt_preproc, wt_preproc_prcp, wt_preproc_prcp_eobs = \
+                load_wt_preprocessors(dataset_name, preproc_variables_dict)
+
+            # calculate lwt
+            lwt = wt_algorithm(wt_preproc, dataset_name)
+
+            era5_ancestors, eobs_ancestors = get_ancestors_era5_eobs(
+                dataset_name, preproc_variables_dict)
+
+            # calculate simplified lwt based on precipitation patterns
+            _ = calc_slwt_obs(
+                    cfg,
+                    lwt,
+                    wt_preproc_prcp,
+                    dataset_name,
+                    era5_ancestors,
+                )
+            _ = calc_slwt_obs(
+                cfg,
+                lwt,
+                wt_preproc_prcp_eobs,
+                'E-OBS',
+                eobs_ancestors,
+            )
+
+            # write only lwt to file
+            write_lwt_to_file(cfg, lwt, wt_preproc, dataset_name)
+
+            # load wt files
+            wt_cubes = load_wt_files(f'{work_dir}/{dataset_name}.nc',
+                                     only_lwt=True)
+
+            var_dict = get_looping_dict(
+                dataset_vars
+            )  # dataset_vars is list of variables for dataset dataset_name
+
+            if plotting:
+                # plot means
+                for var_name, var_data in var_dict.items():
+                    data_info = {'dataset': dataset_name,
+                                 'var': var_name,
+                                 'preproc_path': var_data[1]}
+                    plot_means(cfg, var_data[0], wt_cubes, data_info,
+                               only_lwt=True)
+                plot_seasonal_occurrence(cfg, wt_cubes, dataset_name)
+        else:
+            if dataset_name == 'E-OBS':
+                continue
+            wt_preproc = iris.load_cube(dataset_vars[0].get('filename'))
+
+            output_file_path = get_model_output_filepath(
+                dataset_name, dataset_vars)[1]
+
+            # calculate weathertypes
+            calc_lwt_model(cfg, wt_preproc, dataset_name, output_file_path)
+
+            # load wt files
+            wt_cubes = load_wt_files(f'{work_dir}/{dataset_name}.nc',
+                                     only_lwt=True)
+
+            var_dict = get_looping_dict(
+                dataset_vars
+            )  # dataset_vars is list of variables for dataset dataset_name
+
+            if plotting:
+                # plot means
+                ensemble = dataset_vars[0].get('ensemble')
+                for var_name, var_data in var_dict.items():
+                    data_info = {'dataset': dataset_name,
+                                 'var': var_name,
+                                 'preproc_path': var_data[1],
+                                 'ensemble': ensemble}
+                    plot_means(cfg, var_data[0], wt_cubes, data_info,
+                               only_lwt=True)
+                plot_seasonal_occurrence(cfg, wt_cubes, dataset_name,
+                                         ensemble=ensemble)
+
+
+def run_my_diagnostic(cfg: dict):
+    '''
+    Arguments:
+        cfg - nested dictionary of metadata
+
+    Returns:
+        string; runs the user diagnostic
+
+    '''
+    # assemble the data dictionary keyed by dataset name
+    # this makes use of the handy group_metadata function that
+    # orders the data by 'dataset'; the resulting dictionary is
+    # keyed on datasets e.g. dict = {'MPI-ESM-LR': [var1, var2...]}
+    # where var1, var2 are dicts holding all needed information per variable
+    automatic_slwt = cfg.get('automatic_slwt')
+
+    if automatic_slwt:
+        run_automatic_slwt(cfg)
+    # if automatic_slwt is false, and predefined_slwt is false,
+    # just write selected pairs to file
+    elif not automatic_slwt:
+        run_lwt(cfg)
+
+
+if __name__ == '__main__':
+
+    # always use run_diagnostic() to get the config (the preprocessor
+    # nested dictionary holding all the needed information)
+    with run_diagnostic() as config:
+        # list here the functions that need to run
+        run_my_diagnostic(config)