diff --git a/.gitignore b/.gitignore index 7c783dce5a..abec495ffb 100644 --- a/.gitignore +++ b/.gitignore @@ -112,3 +112,7 @@ esmvaltool/utils/testing/regression/.service/ # Temporary imagehash files imagehashes_*.yml + + +#Ignore vscode AI rules +.github/instructions/codacy.instructions.md diff --git a/doc/sphinx/source/recipes/figures/weathertyping/ERA5__lwt_rel_occurrence_1958-2014.png b/doc/sphinx/source/recipes/figures/weathertyping/ERA5__lwt_rel_occurrence_1958-2014.png new file mode 100755 index 0000000000..d42fbf8784 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/weathertyping/ERA5__lwt_rel_occurrence_1958-2014.png differ diff --git a/doc/sphinx/source/recipes/figures/weathertyping/correlation_matrix_E-OBS_1958-2014.png b/doc/sphinx/source/recipes/figures/weathertyping/correlation_matrix_E-OBS_1958-2014.png new file mode 100755 index 0000000000..6288613b4a Binary files /dev/null and b/doc/sphinx/source/recipes/figures/weathertyping/correlation_matrix_E-OBS_1958-2014.png differ diff --git a/doc/sphinx/source/recipes/figures/weathertyping/lwt_1_ERA5__psl_mean_1958-2014.png b/doc/sphinx/source/recipes/figures/weathertyping/lwt_1_ERA5__psl_mean_1958-2014.png new file mode 100755 index 0000000000..b2cba03df6 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/weathertyping/lwt_1_ERA5__psl_mean_1958-2014.png differ diff --git a/doc/sphinx/source/recipes/figures/weathertyping/lwt_1_TaiESM1_r1i1p1f1_psl_mean_1950-2014.png b/doc/sphinx/source/recipes/figures/weathertyping/lwt_1_TaiESM1_r1i1p1f1_psl_mean_1950-2014.png new file mode 100755 index 0000000000..c8f8ad5fda Binary files /dev/null and b/doc/sphinx/source/recipes/figures/weathertyping/lwt_1_TaiESM1_r1i1p1f1_psl_mean_1950-2014.png differ diff --git a/doc/sphinx/source/recipes/figures/weathertyping/slwt_EOBS_4_ERA5__psl_mean_1958-2014.png b/doc/sphinx/source/recipes/figures/weathertyping/slwt_EOBS_4_ERA5__psl_mean_1958-2014.png new file mode 100755 index 0000000000..1f15102db7 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/weathertyping/slwt_EOBS_4_ERA5__psl_mean_1958-2014.png differ diff --git a/doc/sphinx/source/recipes/index.rst b/doc/sphinx/source/recipes/index.rst index 77a9815107..2a6afd8124 100644 --- a/doc/sphinx/source/recipes/index.rst +++ b/doc/sphinx/source/recipes/index.rst @@ -67,6 +67,7 @@ Atmosphere recipe_radiation_budget recipe_aod_aeronet_assess recipe_surface_trace_gas + recipe_weathertyping Climate metrics ^^^^^^^^^^^^^^^ diff --git a/doc/sphinx/source/recipes/recipe_weathertyping.rst b/doc/sphinx/source/recipes/recipe_weathertyping.rst new file mode 100755 index 0000000000..f7fab453c2 --- /dev/null +++ b/doc/sphinx/source/recipes/recipe_weathertyping.rst @@ -0,0 +1,116 @@ +.. _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 area 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, pr + + +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.* + +This recipe currently only works with the following reanalysis and observation datasets: + +* E-OBS: European Climate Assessment & Dataset gridded daily precipitation sum +* ERA5: ECMWF reanalysis + +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_1958-2014.png + :align: center + + PSL mean map of Lamb WT 1 for ERA5. + +.. _fig_weathertyping_2: +.. figure:: /recipes/figures/weathertyping/lwt_1_TaiESM1_r1i1p1f1_psl_mean_1950-2014.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_1958-2014.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_1958-2014.png + :align: center + + Heatmap of correlation values for Lamb WTs 1-27. + +.. _fig_weathertyping_5: +.. figure:: /recipes/figures/weathertyping/ERA5__lwt_rel_occurrence_1958-2014.png + :align: center + + Stackplot of seasonal relative occurrences of each WT for ERA5. diff --git a/esmvaltool/config-references.yml b/esmvaltool/config-references.yml index cda65ff656..939696123f 100644 --- a/esmvaltool/config-references.yml +++ b/esmvaltool/config-references.yml @@ -322,6 +322,10 @@ authors: name: Juckes, Martin institute: BADC, UK orcid: + jury_martin: + name: Jury, Martin + institute: WEGC, Austria + orcid: https://orcid.org/0000-0003-0590-7843 kadygrov_nikolay: name: Kadygrov, Nikolay institute: IPSL, France @@ -355,6 +359,11 @@ authors: name: Krasting, John institute: NOAA, USA orcid: https://orcid.org/0000-0002-4650-9844 + kroissenbrunner_thomas: + name: Kroissenbrunner, Thomas + institute: WEGC, Austria + github: thomaskroi1996 + orcid: kuehbacher_birgit: name: Kuehbacher, Birgit institute: DLR, Germany @@ -881,6 +890,7 @@ projects: ipcc_ar6: IPCC AR6 WG1 contributions isenes3: EU H2020 Infrastructure for the European Network for Earth System Modelling - Phase 3 primavera: EU H2020 project PRIMAVERA + preval: PREVAL ÖKS qa4ecv: QA4ECV russell_project: US Clivar, Ocean Carbon Biogeochemistry, ESGF, NOAA, NSF, NASA, US Antarctic program trr181: DFG Project TRR-181, Energy transfers in Atmosphere and Ocean @@ -888,6 +898,7 @@ projects: usmile: ERC Synergy Grant USMILE + realms: aerosol: aerosol atmos: atmosphere diff --git a/esmvaltool/diag_scripts/weathertyping/__init__.py b/esmvaltool/diag_scripts/weathertyping/__init__.py new file mode 100644 index 0000000000..e31afa296a --- /dev/null +++ b/esmvaltool/diag_scripts/weathertyping/__init__.py @@ -0,0 +1 @@ +"""Initialize the ESMValTool weathertyping package.""" diff --git a/esmvaltool/diag_scripts/weathertyping/calc_utils.py b/esmvaltool/diag_scripts/weathertyping/calc_utils.py new file mode 100644 index 0000000000..86e8ba2584 --- /dev/null +++ b/esmvaltool/diag_scripts/weathertyping/calc_utils.py @@ -0,0 +1,1104 @@ +"""Utility functions for calculations.""" + +import json +import logging +import warnings +from pathlib import Path +from typing import Final + +import dask.array as da +import iris +import iris.analysis.cartography +import numpy as np +import pandas as pd +from plot_utils import plot_corr_rmse_heatmaps, plot_maps +from wt_utils import ( + check_mapping_dict_format, + get_driver, + get_mapping_dict, + get_provenance_record, + log_provenance, + map_lwt_to_slwt, + reverse_convert_dict, + write_corr_rmse_to_csv, + write_mapping_dict, +) + +iris.FUTURE.datum_support = True + +logger = logging.getLogger(Path(__file__).name) + +# Ignoring a warning that is produced when selecting timesteps of a weathertype +warnings.filterwarnings("ignore", ".*Collapsing a non-contiguous coordinate*") + + +# defining constants for direction calculations +DIR_0: Final = 0.0 +DIR_22_5: Final = 22.5 +DIR_67_5: Final = 67.5 +DIR_112_5: Final = 112.5 +DIR_157_5: Final = 157.5 +DIR_202_5: Final = 202.5 +DIR_247_5: Final = 247.5 +DIR_292_5: Final = 292.5 +DIR_337_5: Final = 337.5 + +# defining constant for flow limit +FLOW_LIMIT: Final = 6 + + +def calc_slwt_obs( + cfg: dict, + lwt: np.array, + cube: iris.cube.Cube, + dataset: str, + ancestors: list, + timerange: str, +) -> np.array: + """Calculate simplified weathertypes for observational data. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + lwt + Array of Lamb WT + cube + Cube of psl data + dataset + Name of dataset + ancestors + List of ancestors + timerange + Time range for the calculation + + Returns + ------- + np.array: Simplified Lamb Weathertypes + + """ + logger.info("Calculating simplified Lamb Weathertypes for %s", dataset) + + wt_data_prcp = [] + for wt_ in range(1, 28): + target_indices = da.where(lwt == wt_)[0].compute() + if target_indices.size == 0: + logger.info( + "calc_slwt_obs - CAUTION: Skipped wt %s \ + for dataset %s!", + wt_, + dataset, + ) + continue + dates = [ + cube.coord("time").units.num2date(cube.coord("time").points[i]) + for i in target_indices + ] + if dataset == "E-OBS": + extracted_cube = cube[target_indices] + else: + extracted_cube = cube.extract( + iris.Constraint(time=lambda t, d=dates: t.point in d), + ) + wt_cube_mean = extracted_cube.collapsed("time", iris.analysis.MEAN) + + data = wt_cube_mean.data + + if isinstance(data, np.ma.MaskedArray): + wt_data_prcp.append(data.compressed()) + else: + # if not MaskedArray, flatten should be fine + wt_data_prcp.append(data.flatten()) + + selected_pairs = process_prcp_mean( + cfg, + wt_data_prcp, + dataset, + timerange, + ) + + with Path( + f"{cfg.get('work_dir')}/wt_selected_pairs_{dataset}.json", + ).open("w", encoding="utf-8") as file: + json.dump(selected_pairs, file) + + mapping_dict = get_mapping_dict(selected_pairs) + + write_mapping_dict(cfg.get("work_dir"), dataset, mapping_dict) + + provenance_record = get_provenance_record( + "Lamb Weathertypes", + ancestors, + ["Lamb Weathertypes"], + ) + + log_provenance( + f"{cfg.get('work_dir')}/wt_selected_pairs_{dataset}", + cfg, + provenance_record, + ) + + return map_lwt_to_slwt(lwt, mapping_dict) + + +def calc_const(): + """Calculate constants for weathertyping algorithm. + + Eq. taken from: 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 + + Returns + ------- + tuple + The four constants needed for WT calculation. + """ + const1 = 1 / np.cos(45 * np.pi / 180) + const2 = np.sin(45 * np.pi / 180) / np.sin(40 * np.pi / 180) + const3 = np.sin(45 * np.pi / 180) / np.sin(50 * np.pi / 180) + const4 = 1 / (2 * np.cos(45 * np.pi / 180) ** 2) + + return const1, const2, const3, const4 + + +def calc_westerly_flow(cube: iris.cube.Cube) -> np.array: + """Calculate the westerly flow over area. + + Eq. taken from: 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 + + Parameters + ---------- + cube + Cube of psl data. + + Returns + ------- + np.array + westerly flow + """ + return 1 / 2 * (cube.data[:, 1, 2] + cube.data[:, 1, 4]) - 1 / 2 * ( + cube.data[:, 3, 2] + cube.data[:, 3, 4] + ) + + +def calc_southerly_flow(cube: iris.cube.Cube, const1: float) -> np.array: + """Calculate the southerly flow over area. + + Eq. taken from: 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 + + Parameters + ---------- + cube + cube of psl data + const1 + const1 + + Returns + ------- + np.array + southerly flow + """ + return const1 * ( + 1 + / 4 + * (cube.data[:, 3, 4] + 2 * cube.data[:, 2, 4] + cube.data[:, 1, 4]) + - 1 + / 4 + * (cube.data[:, 3, 2] + 2 * cube.data[:, 2, 2] + cube.data[:, 1, 2]) + ) + + +def calc_resultant_flow( + westerly_flow: np.array, + southerly_flow: np.array, +) -> np.array: + """Calculate the resultant flow. + + Eq. taken from: 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 + + Parameters + ---------- + westerly_flow + westerly flow + southerly_flow + southerly flow + + Returns + ------- + np.array + resultant flow + """ + return (southerly_flow**2 + westerly_flow**2) ** (1 / 2) + + +def calc_westerly_shear_velocity( + cube: iris.cube.Cube, + const2: float, + const3: float, +) -> np.array: + """Calculate westerly shear velocity. + + Eq. taken from: 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 + + Parameters + ---------- + cube + cube of psl data + const2 + const2 + const3 + const3 + + Returns + ------- + np.array + westerly shear velocity + """ + return const2 * ( + 1 / 2 * (cube.data[:, 0, 2] + cube.data[:, 0, 4]) + - 1 / 2 * (cube.data[:, 2, 2] + cube.data[:, 2, 4]) + ) - const3 * ( + 1 / 2 * (cube.data[:, 2, 2] + cube.data[:, 2, 4]) + - 1 / 2 * (cube.data[:, 4, 2] + cube.data[:, 4, 4]) + ) + + +def calc_southerly_shear_velocity( + cube: iris.cube.Cube, + const4: float, +) -> np.array: + """Calculate southerly shear velocity. + + Eq. taken from: 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 + + Parameters + ---------- + cube + cube of psl data + const4 + const + + Returns + ------- + np.array + southerly shear velocity + """ + return const4 * ( + 1 + / 4 + * (cube.data[:, 3, 6] + 2 * cube.data[:, 2, 6] + cube.data[:, 1, 6]) + - 1 + / 4 + * (cube.data[:, 3, 4] + 2 * cube.data[:, 2, 4] + cube.data[:, 1, 4]) + - 1 + / 4 + * (cube.data[:, 3, 2] + 2 * cube.data[:, 2, 2] + cube.data[:, 1, 2]) + + 1 + / 4 + * (cube.data[:, 3, 0] + 2 * cube.data[:, 2, 0] + cube.data[:, 1, 0]) + ) + + +def calc_total_shear_velocity( + westerly_shear_velocity: np.array, + southerly_shear_velocity: np.array, +) -> np.array: + """Calculate total shear velocity. + + Eq. taken from: 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 + + Parameters + ---------- + westerly_shear_velocity + westerly shear velocity + southerly_shear_velocity + southerly shear velocity + + Returns + ------- + np.array + total shear velocity + """ + return westerly_shear_velocity + southerly_shear_velocity + + +def lamp_pure_directional_type( + direction: float, + weathertypes: np.array, + i: int, +) -> int: + """Calculate Lamp pure directional weathertype. + + Parameters + ---------- + direction + direction in degrees + + Returns + ------- + int + Lamb weathertype + """ + if direction >= DIR_337_5 or direction < DIR_22_5: + weathertypes[i] = 1 + if DIR_22_5 <= direction < DIR_67_5: + weathertypes[i] = 2 + if DIR_67_5 <= direction < DIR_112_5: + weathertypes[i] = 3 + if DIR_112_5 <= direction < DIR_157_5: + weathertypes[i] = 4 + if DIR_157_5 <= direction < DIR_202_5: + weathertypes[i] = 5 + if DIR_202_5 <= direction < DIR_247_5: + weathertypes[i] = 6 + if DIR_247_5 <= direction < DIR_292_5: + weathertypes[i] = 7 + if DIR_292_5 <= direction < DIR_337_5: + weathertypes[i] = 8 + + +def lamp_synoptic_directional_type( + direction: float, + weathertypes: np.array, + i: int, +) -> int: + """Calculate Lamb synoptic/directional hybrid weathertype. + + Parameters + ---------- + direction + direction in degrees + + Returns + ------- + int + Lamb weathertype + """ + if direction >= DIR_337_5 or direction < DIR_22_5: + weathertypes[i] = 11 + if DIR_22_5 <= direction < DIR_67_5: + weathertypes[i] = 12 + if DIR_67_5 <= direction < DIR_112_5: + weathertypes[i] = 13 + if DIR_112_5 <= direction < DIR_157_5: + weathertypes[i] = 14 + if DIR_157_5 <= direction < DIR_202_5: + weathertypes[i] = 15 + if DIR_202_5 <= direction < DIR_247_5: + weathertypes[i] = 16 + if DIR_247_5 <= direction < DIR_292_5: + weathertypes[i] = 17 + if DIR_292_5 <= direction < DIR_337_5: + weathertypes[i] = 18 + + +def lamb_synoptic_directional_type_zlt0( + direction: float, + weathertypes: np.array, + i: int, +) -> int: + """Calculate Lamb synoptic/directional hybrid weathertype with z <0. + + Parameters + ---------- + direction + direction in degrees + + Returns + ------- + int + Lamb weathertype + """ + if direction >= DIR_337_5 or direction < DIR_22_5: + weathertypes[i] = 19 + if DIR_22_5 <= direction < DIR_67_5: + weathertypes[i] = 20 + if DIR_67_5 <= direction < DIR_112_5: + weathertypes[i] = 21 + if DIR_112_5 <= direction < DIR_157_5: + weathertypes[i] = 22 + if DIR_157_5 <= direction < DIR_202_5: + weathertypes[i] = 23 + if DIR_202_5 <= direction < DIR_247_5: + weathertypes[i] = 24 + if DIR_247_5 <= direction < DIR_292_5: + weathertypes[i] = 25 + if DIR_292_5 <= direction < DIR_337_5: + weathertypes[i] = 26 + + +def wt_algorithm(cube: iris.cube.Cube, dataset: str) -> np.array: + """Algorithm to calculate Lamb weathertypes. + + Eq. taken from: 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 + + Parameters + ---------- + cube + Cube of psl data + dataset + Name of dataset + + Returns + ------- + np.array + Lamb weathertypes + """ + # lats and lons corresponding to datapoints + # 55, 5 -> 1 + # 55, 15 -> 2 + # 50, -5 -> 3 + # 50, 5 -> 4 + # 50, 15 -> 5 + # 50, 25 -> 6 + # 45, -5 -> 7 + # 45, 5 -> 8 + # 45, 15 -> 9 + # 45, 25 -> 10 + # 40, -5 -> 11 + # 40, 5 -> 12 + # 40, 15 -> 13 + # 40, 25 -> 14 + # 35, 5 -> 15 + # 35, 15 -> 16 + + # lons: -5, 0, 5, 10, 15, 20, 25 + # lats: 35, 40, 45, 50, 55 + logger.info("Calculating Lamb Weathertypes for %s", dataset) + + const1, const2, const3, const4 = calc_const() + + # westerly flow + westerly_flow = calc_westerly_flow(cube) + # southerly flow + southerly_flow = calc_southerly_flow(cube, const1) + # resultant flow + total_flow = calc_resultant_flow(westerly_flow, southerly_flow) + # westerly shear vorticity + westerly_shear_velocity = calc_westerly_shear_velocity( + cube, + const2, + const3, + ) + # southerly shear vorticity + southerly_shear_velocity = calc_southerly_shear_velocity(cube, const4) + # total shear vorticity + total_shear_velocity = calc_total_shear_velocity( + westerly_shear_velocity, + southerly_shear_velocity, + ) + + weathertypes = np.zeros(len(total_shear_velocity)) + + for i, z_i in enumerate(total_shear_velocity): + direction = ( + np.arctan(westerly_flow[i] / southerly_flow[i]) * 180 / np.pi + ) # deg + if southerly_flow[i] >= 0: + direction += 180 # deg + + if direction < DIR_0: + direction += 360 # deg + + # Lamb pure directional type + if abs(z_i) < total_flow[i]: + # writes directly to weathertypes array + lamp_pure_directional_type(direction, weathertypes, i) + # Lamb pure cyclonic and anticyclonic type + elif (2 * total_flow[i]) < abs(z_i): + if z_i > 0: + weathertypes[i] = 9 + + elif z_i < 0: + weathertypes[i] = 10 + # Lamb synoptic/direction hybrid types + elif total_flow[i] < abs(z_i) < (2 * total_flow[i]): + if z_i > 0: + # writes directly to weathertypes array + lamp_synoptic_directional_type(direction, weathertypes, i) + + elif z_i < 0: + # writes directly to weathertypes array + lamb_synoptic_directional_type_zlt0(direction, weathertypes, i) + # light indeterminate flow, corresponding to Lamb unclassified type U + elif abs(z_i) < FLOW_LIMIT and total_flow[i] < FLOW_LIMIT: + weathertypes[i] = 27 + + return weathertypes + + +def calc_lwt_slwt_model( + cfg: dict, + cube: iris.cube.Cube, + data_info: dict, + predefined_slwt: dict | None = None, +): + """Calculate Lamb as well as simplified weathertypes for model and write them to file. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + cube + PSL field of dataset + data_info + Dictionary with info to dataset + predefined_slwt + If None, automatic_slwt will be used. + If dict, this mapping dict will be used. + (see recipe option predefined_slwt) + """ + driver = get_driver(data_info) + + if not Path( + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}", + ).exists(): + Path( + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}", + ).mkdir( + parents=True, + exist_ok=True, + ) + + lwt = wt_algorithm(cube, data_info.get("dataset")) + + time_points = cube.coord("time").units.num2date(cube.coord("time").points) + + logger.info( + "Calculating simplified Lamb Weathertypes for %s %s", + data_info.get("dataset"), + data_info.get("ensemble", ""), + ) + + if not predefined_slwt: + with Path( + f"{cfg.get('work_dir')}/wt_mapping_dict_ERA5.json", + ).open("r", encoding="utf-8") as file: + mapping_dict_era5_f = json.load(file) + mapping_dict_era5 = reverse_convert_dict(mapping_dict_era5_f) + + with Path( + f"{cfg.get('work_dir')}/wt_mapping_dict_E-OBS.json", + ).open("r", encoding="utf-8") as file: + mapping_dict_eobs_f = json.load(file) + mapping_dict_eobs = reverse_convert_dict(mapping_dict_eobs_f) + + slwt_era5 = map_lwt_to_slwt(lwt, mapping_dict_era5) + slwt_eobs = map_lwt_to_slwt(lwt, mapping_dict_eobs) + + print("ERA5: ", slwt_era5) + print("EOBS: ", slwt_eobs) + else: + predefined_slwt = check_mapping_dict_format(predefined_slwt) + write_mapping_dict(cfg.get("work_dir"), "ERA5", predefined_slwt) + write_mapping_dict(cfg.get("work_dir"), "E-OBS", predefined_slwt) + slwt_era5 = map_lwt_to_slwt(lwt, predefined_slwt) + slwt_eobs = map_lwt_to_slwt(lwt, predefined_slwt) + + wt_cube = iris.cube.CubeList() + wt_cube.append(iris.cube.Cube(lwt, long_name="lwt")) + wt_cube.append(iris.cube.Cube(slwt_era5, long_name="slwt_era5")) + wt_cube.append(iris.cube.Cube(slwt_eobs, long_name="slwt_eobs")) + + wt_cube[0].add_dim_coord(cube.coord("time"), 0) + wt_cube[1].add_dim_coord(cube.coord("time"), 0) + wt_cube[2].add_dim_coord(cube.coord("time"), 0) + + iris.save( + wt_cube, + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}/{data_info.get('dataset')}{driver}_" + f"{data_info.get('ensemble', '')}_{data_info.get('timerange')}.nc", + ) + + # write to csv file + d = { + "date": time_points[:], + "lwt": np.int8(lwt), + "slwt_ERA5": np.int8(slwt_era5), + "slwt_EOBS": np.int8(slwt_eobs), + } + df = pd.DataFrame(data=d) + df.to_csv( + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}/{data_info.get('dataset')}{driver}_{data_info.get('ensemble', '')}_" + f"{data_info.get('timerange')}.csv", + index=False, + ) + + ancestors = [ + data_info.get("preproc_path"), + f"{cfg.get('work_dir')}/wt_mapping_dict_ERA5.json", + f"{cfg.get('work_dir')}/wt_mapping_dict_E-OBS.json", + ] + provenance_record = get_provenance_record( + "Lamb Weathertypes", + ancestors, + ["Lamb Weathertypes"], + ) + + log_provenance( + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}/{data_info.get('dataset')}_{data_info.get('ensemble', '')}", + cfg, + provenance_record, + ) + + +def rmse(subarray1: np.array, subarray2: np.array) -> np.array: + """Calculate root mean square error between two arrays.""" + return np.sqrt(np.mean((subarray1 - subarray2) ** 2)) + + +def process_prcp_mean( + cfg: dict, + data: np.array, + dataset: str, + timerange: str, +) -> list: + """Return which weathertypes can be grouped together for a certain precipitation pattern. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + data + Array of precipitation means for each WT + dataset + Name of dataset + timerange + Time range of dataset + + Returns + ------- + list + Selected pairs of WT. This is passed to get_mapping_dict + """ + logger.info("Calculating corr and rsme matrices for %s", dataset) + + selected_pairs = [] + pattern_correlation_matrix = np.ma.corrcoef(data) + n = len(data) + rmse_matrix = np.zeros((n, n)) + + for i in range(n): + for j in range(i + 1, n): + rmse_matrix[i][j] = rmse(data[i], data[j]) + rmse_matrix[j][i] = rmse_matrix[i][j] + if pattern_correlation_matrix[i][j] >= cfg.get( + "correlation_threshold", + ) and rmse_matrix[i][j] <= cfg.get("rmse_threshold"): + selected_pairs.append( + (i + 1, j + 1), + ) + + # write matrices to csv + write_corr_rmse_to_csv( + cfg, + pattern_correlation_matrix, + rmse_matrix, + dataset, + ) + # plot heatmaps for matrices + plot_corr_rmse_heatmaps( + cfg, + pattern_correlation_matrix, + rmse_matrix, + dataset, + timerange, + ) + + return selected_pairs + + +def calc_wt_means( + cfg: dict, + cube: iris.cube.Cube, + wt_cubes: iris.cube.CubeList, + data_info: dict, +): + """Calculate means for psl, tas or pr for weathertypes. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + cube + Cube with variable data + wt_cubes + List of cubes of lwt, slwt_ERA5 and slwt_EOBS + data_info + Dictionary with info to dataset + """ + var_name = data_info.get("var") + wt_string = data_info.get("wt_string") + driver = get_driver(data_info) + + logger.info( + "Calculating %s %s means for %s", + data_info.get("dataset"), + var_name, + wt_string, + ) + + wt_array, tcoord = get_wt_array(wt_string, wt_cubes) + counts = da.bincount(wt_array.flatten().astype(int)).compute() + + for wt in range(1, len(counts)): + target_indices = da.where(wt_array == wt)[0].compute() + if target_indices.size == 0: + logger.info( + "calc_wt_means - CAUTION: Skipped %s %s \ + for dataset %s!", + wt_string, + wt, + data_info.get("dataset"), + ) + continue + dates = [ + tcoord.units.num2date(tcoord.points[i]) for i in target_indices + ] + extracted_cube = cube.extract( + iris.Constraint(time=lambda t, d=dates: t.point in d), + ) + wt_cube_mean = extracted_cube.collapsed("time", iris.analysis.MEAN) + plot_maps(wt, cfg, wt_cube_mean, data_info, "mean") + + ancestors = [ + f"{data_info.get('preproc_path')}", + f"{cfg.get('work_dir')}/ERA5.nc", + ] + provenance_record = get_provenance_record( + f"{var_name} means for {wt_string}, wt: {wt}, ", + ancestors, + [var_name], + ["map"], + ["mean"], + ) + + local_path = f"{cfg.get('plot_dir')}/mean" + + log_provenance( + f"{local_path}/{wt_string}_{wt}{driver}_" + f"{data_info.get('dataset')}_{data_info.get('ensemble', '')}" + f"_{var_name}_mean_{data_info.get('timerange')}", + cfg, + provenance_record, + ) + + +def get_wt_array(wt_string: str, wt_cubes: iris.cube.CubeList) -> tuple: + """Get weathertype array and time coordinate based on wt_string. + + Parameters + ---------- + wt_string + string for weathertype selection + wt_cubes + list of weathertype cubes + + Raises + ------ + NameError + if wt_array does not exist for the given wt_string + + Returns + ------- + tuple + weathertype array and time coordinate + """ + if wt_string == "slwt_ERA5": + slwt_era5_cube = wt_cubes[1] + tcoord = slwt_era5_cube.coord("time") + wt_array = slwt_era5_cube.core_data()[:] + elif wt_string == "slwt_EOBS": + slwt_eobs_cube = wt_cubes[2] + tcoord = slwt_eobs_cube.coord("time") + wt_array = slwt_eobs_cube.core_data()[:] + elif wt_string == "lwt": + lwt_cube = wt_cubes[0] + tcoord = lwt_cube.coord("time") + wt_array = lwt_cube.core_data()[:] + else: + error_str = "wt_array does not exist for calc_wt_means, calc_wt_anomalies or calc_wt_std." + raise NameError( + error_str, + ) + + return wt_array, tcoord + + +def calc_wt_anomalies( + cfg: dict, + cube: iris.cube.Cube, + wt_cubes: iris.cube.CubeList, + data_info: dict, +): + """Calculate anomalies for psl, tas and pr for weathertypes. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + cube + Cube with variable data + wt_cubes + List of cubes of lwt, slwt_ERA5 and slwt_EOBS + data_info + Dictionary with info to dataset + """ + var_name = data_info.get("var_name") + wt_string = data_info.get("wt_string") + + logger.info( + "Calculating %s %s anomalies for %s", + data_info.get("dataset"), + var_name, + wt_string, + ) + + wt_array, tcoord = get_wt_array(wt_string, wt_cubes) + + for wt in range(1, len(np.unique(wt_array)) + 1): + target_indices = da.where(wt_array == wt).compute() + if target_indices.size == 0: + logger.info( + "calc_wt_anomalies - CAUTION: Skipped wt %s \ + for dataset %s!", + wt, + data_info.get("dataset"), + ) + continue + dates = [ + tcoord.units.num2date(tcoord.points[i]) for i in target_indices + ] + extracted_cube = cube.extract( + iris.Constraint(time=lambda t, d=dates: t.point in d[0]), + ) + wt_cube_mean = extracted_cube.collapsed("time", iris.analysis.MEAN) + plot_maps( + wt, + cfg, + cube.collapsed("time", iris.analysis.MEAN) - wt_cube_mean, + data_info, + "anomaly", + ) + + ancestors = [ + f"{data_info.get('preproc_path')}", + f"{cfg.get('work_dir')}/ERA5.nc", + ] + provenance_record = get_provenance_record( + f"{var_name} anomaly for, wt: {wt} \ + {wt_string}", + ancestors, + [var_name], + ["map"], + ["anomaly"], + ) + + log_provenance( + f"{cfg.get('plot_dir')}/anomaly/{wt_string}_{wt}_{data_info.get('dataset')}_{data_info.get('ensemble', '')}" + f"_{var_name}_anomaly__{data_info.get('timerange')}", + cfg, + provenance_record, + ) + + +def calc_wt_std( + cfg: dict, + cube: iris.cube.Cube, + wt_cubes: iris.cube.CubeList, + data_info: dict, +): + """Calculate standard deviation for psl, tas and pr for weathertypes. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + cube + Cube with variable data + wt_cubes + List of cubes of lwt, slwt_ERA5 and slwt_EOBS + data_info + Dictionary with info to dataset + """ + var_name = data_info.get("var_name") + wt_string = data_info.get("wt_string") + + logger.info( + "Calculating %s %s standard deviation for %s", + data_info.get("dataset"), + var_name, + wt_string, + ) + + wt_array, tcoord = get_wt_array(wt_string, wt_cubes) + + for wt in range(1, len(np.unique(wt_array)) + 1): + target_indices = da.where(wt_array == wt).compute() + if target_indices.size == 0: + logger.info( + "calc_slwt_obs - CAUTION: Skipped wt %s \ + for dataset %s!", + wt, + data_info.get("dataset"), + ) + continue + dates = [ + tcoord.units.num2date(tcoord.points[i]) for i in target_indices + ] + extracted_cube = cube.extract( + iris.Constraint(time=lambda t, d=dates: t.point in d[0]), + ) + wt_cube_std = extracted_cube.collapsed("time", iris.analysis.STD_DEV) + plot_maps(wt, cfg, wt_cube_std, data_info, "stddev") + + ancestors = [ + f"{data_info.get('preproc_path')}", + f"{cfg.get('work_dir')}/ERA5.nc", + ] + provenance_record = get_provenance_record( + f"{var_name} standard, wt: {wt}\ + deviation for \ + {wt_string}", + ancestors, + [var_name], + ["map"], + ["stddev"], + ) + + local_path = f"{cfg.get('plot_dir')}/stddev" + + log_provenance( + f"{local_path}/{wt_string}_{wt}_{data_info.get('dataset')}_{data_info.get('ensemble', '')}" + f"_{var_name}_stddev_{data_info.get('timerange')}", + cfg, + provenance_record, + ) + + +def calc_lwt_model(cfg: dict, cube: iris.cube.Cube, data_info: dict): + """Calculate Lamb weathertypes for models. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + cube + Cube with variable data + data_info + Dictionary with info to dataset + """ + driver = get_driver(data_info) + + if not Path( + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}", + ).exists(): + Path( + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}", + ).mkdir(parents=True, exist_ok=True) + + wt = wt_algorithm(cube, data_info.get("dataset")) + + time_points = cube.coord("time").units.num2date(cube.coord("time").points) + + wt_cube = iris.cube.CubeList() + wt_cube.append(iris.cube.Cube(wt, long_name="lwt")) + + logger.info( + "Writing Lamb Weathertype for %s \ + to file %s.nc", + data_info.get("dataset"), + data_info.get("dataset"), + ) + + wt_cube[0].add_dim_coord(cube.coord("time"), 0) + + iris.save( + wt_cube, + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}/{data_info.get('dataset')}{driver}" + f"_{data_info.get('ensemble', '')}_{data_info.get('timerange')}.nc", + ) + + # write to csv file + d = {"date": time_points[:], "lwt": np.int8(wt)} + df = pd.DataFrame(data=d) + df.to_csv( + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}/{data_info.get('dataset')}{driver}_{data_info.get('ensemble', '')}_" + f"{data_info.get('timerange')}.csv", + index=False, + ) + + ancestors = [ + f"{cfg.get('work_dir')}/wt_mapping_dict_ERA5.json", + f"{cfg.get('work_dir')}/wt_mapping_dict_E-OBS.json", + ] + provenance_record = get_provenance_record( + "Lamb Weathertypes", + ancestors, + ["Lamb Weathertypes"], + ) + + log_provenance( + f"{cfg.get('work_dir')}/{data_info.get('output_file_path')}/{data_info.get('dataset')}_{data_info.get('ensemble', '')}", + cfg, + provenance_record, + ) + + +def plot_means( + cfg: dict, + preproc_var: np.array, + wt_cubes: iris.cube.Cube, + data_info: dict, + mode: str = "slwt", +): + """Plot means, anomalies and standard deviations. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + preproc_var + Preprocessed variable cube + wt_cubes + List of cubes of lwt, slwt_ERA5 and slwt_EOBS + data_info + Dictionary with info to dataset + mode + Mode of weathertype calculation, either "slwt" or "lwt". Defaults to "slwt". + """ + if mode == "slwt": + data_info["wt_string"] = "lwt" + calc_wt_means(cfg, preproc_var, wt_cubes, data_info) + data_info["wt_string"] = "slwt_ERA5" + calc_wt_means(cfg, preproc_var, wt_cubes, data_info) + data_info["wt_string"] = "slwt_EOBS" + calc_wt_means(cfg, preproc_var, wt_cubes, data_info) + elif mode == "lwt": + data_info["wt_string"] = "lwt" + calc_wt_means(cfg, preproc_var, wt_cubes, data_info) + else: + e = "mode must be either 'slwt' or 'lwt'" + raise ValueError(e) diff --git a/esmvaltool/diag_scripts/weathertyping/plot_utils.py b/esmvaltool/diag_scripts/weathertyping/plot_utils.py new file mode 100644 index 0000000000..6f669dfb3f --- /dev/null +++ b/esmvaltool/diag_scripts/weathertyping/plot_utils.py @@ -0,0 +1,556 @@ +"""Utility functions for plotting.""" + +import logging +import warnings +from pathlib import Path + +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import iris +import iris.analysis.cartography +import iris.plot as iplt +import matplotlib.pyplot as plt +import matplotlib.ticker as mticker +import numpy as np +import seaborn as sns +from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER +from matplotlib.colors import ListedColormap +from wt_utils import get_driver, get_provenance_record, log_provenance + +iris.FUTURE.datum_support = True + +logger = logging.getLogger(Path(__file__).name) + +# Ignoring a warning that is produced when selecting timesteps of a weathertype +warnings.filterwarnings("ignore", ".*Collapsing a non-contiguous coordinate*") + + +def generate_grayscale_hex_values(x): + """Generate grayscale values for plotting seasonal occurrences. + + Parameters + ---------- + x + Amount of weathertypes. + + Returns + ------- + list + List of grayscale hex values + """ + grayscale_values = np.linspace(0, 1, x) + + return [ + f"#{int(value * 255):02x}{int(value * 255):02x}{int(value * 255):02x}" + for value in grayscale_values + ] + + +def plot_seasonal_occurrence( + cfg: dict, + wt_cubes: iris.cube.Cube, + data_info: dict, + cube_path: str, +): + """Plot seasonal occurrences of weathertypes. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + wt_cubes + List of cubes of lwt, slwt_ERA5 and slwt_EOBS + data_info + Dictionary with info to dataset + cube_path + Paths to weathertype cubes (ancestors) + """ + driver = get_driver(data_info) + + output_path = f"{cfg['plot_dir']}/seasonal_occurrence" + + if not Path(output_path).exists(): + Path(output_path).mkdir(parents=True, exist_ok=True) + + month_list = [ + "Jan", + "Feb", + "Mar", + "Apr", + "May", + "Jun", + "Jul", + "Aug", + "Sep", + "Oct", + "Nov", + "Dec", + ] + + relative_occurrences = {} + # {wt_string: {month: {wt1: rel_occurrence, wt2: rel_occurrence, ....}}} + # first do absolute occurrence, then relative occurrence + for cube in wt_cubes: + month_dict = {} + for month in range(1, 13): + array = cube.extract( + iris.Constraint(time=iris.time.PartialDateTime(month=month)), + ).data + unique, counts = np.unique(array, return_counts=True) + month_dict[month] = dict( + zip(unique, counts / sum(counts), strict=False), + ) + relative_occurrences[cube.long_name] = month_dict + + for wt_string in relative_occurrences: + wt_numbers = max( + len(value) + for value in relative_occurrences.get(wt_string).values() + ) + wt_stack = np.zeros((wt_numbers, 12)) + for month, month_value in relative_occurrences.get(wt_string).items(): + for wt in month_value: + wt_stack[np.int8(wt - 1), month - 1] = month_value.get(wt) + + y = np.vstack(list(wt_stack)) + + # plot + _, ax_ = plt.subplots(figsize=(10, 10)) + + ax_.set_title( + f"Seasonal occurrence of {wt_string}, \ + {data_info.get('dataset')}, {data_info.get('timerange')}", + ) + + ax_.stackplot( + month_list, + y, + colors=generate_grayscale_hex_values(wt_numbers), + ) + + ax_.legend( + loc="upper center", + bbox_to_anchor=(0.5, -0.05), + fancybox=True, + shadow=True, + ncol=9, + labels=tuple(f"WT {i + 1}" for i in range(wt_numbers)), + ) + + ax_.set( + xlim=(0, 11), + xticks=np.arange(0, 12), + ylim=(0, 1), + yticks=np.arange(0, 1.1, 0.1), + ) + ax_.set_xlabel("Month") + ax_.set_ylabel("Cumulative Relative occurrence") + + plt.savefig( + f"{output_path}/{driver}{data_info.get('dataset')}_{data_info.get('ensemble', '')}_" + f"{wt_string}_rel_occurrence_{data_info.get('timerange')}.png", + ) + plt.savefig( + f"{output_path}/{driver}{data_info.get('dataset')}_{data_info.get('ensemble', '')}_" + f"{wt_string}_rel_occurrence_{data_info.get('timerange')}.pdf", + ) + plt.close() + # ancestors here are just the wt_cubes i guess + ancestors = [ + f"{cube_path}", + ] + provenance_record = get_provenance_record( + f"Seasonal occurrences for {wt_string}, ", + ancestors, + ["wt occurrences"], + ["seas"], + ) + + local_path = f"{cfg.get('plot_dir')}/mean" + + log_provenance( + f"{local_path}/{wt_string}_{driver}_" + f"{data_info.get('dataset')}_{data_info.get('ensemble', '')}" + f"_seasonal_occurrence_{data_info.get('timerange')}", + cfg, + provenance_record, + ) + + +def set_gridlines(ax: plt.Axes): + """Set gridlines for plotting maps. + + Parameters + ---------- + ax + Axes object to draw gridlines on. + + Returns + ------- + gl : (ax.gridlines) + Gridlines object + """ + gl = ax.gridlines( + crs=ccrs.PlateCarree(), + draw_labels=True, + linewidth=0.5, + color="gray", + alpha=0.5, + linestyle="--", + ) + gl.left_labels = True + gl.bottom_labels = True + gl.top_labels = False + gl.right_labels = False + gl.xlines = True + gl.ylocator = mticker.FixedLocator(np.arange(20, 70, 5)) + gl.xlocator = mticker.FixedLocator([-10, -5, 0, 5, 10, 15]) + gl.xformatter = LONGITUDE_FORMATTER + gl.yformatter = LATITUDE_FORMATTER + gl.xlabel_style = {"size": 8, "color": "black"} + gl.ylabel_style = {"color": "black", "size": 8} + + return gl + + +def prepare_plot_title( + data_info: dict, + wt: int, + var_name: str, + mode: str, +) -> str: + """Return formatted plot title. + + Parameters + ---------- + data_info + dict containing info on dataset + wt + weathertype + var_name + name of variable + mode + statistic used + + Returns + ------- + str + title of plot + """ + ensemble = data_info.get("ensemble", "") + timerange = data_info.get("timerange") + dataset = data_info.get("dataset") + if var_name == "pr" and dataset == "ERA5": + title = f"{dataset} {ensemble}, total {var_name} {mode}\n{timerange}, wt: {wt}" + elif var_name == "pr": + title = f"{dataset} {ensemble}, {var_name} flux {mode}\n{timerange}, wt: {wt}" + elif var_name == "tas": + title = f"{dataset} {ensemble}, 1000 hPa {var_name} {mode}\n{timerange}, wt: {wt}" + else: # psl or others + title = ( + f"{dataset} {ensemble}, {var_name} {mode}\n{timerange}, wt: {wt}" + ) + return title + + +def get_unit(var_name: str, dataset: str) -> str: + """Get unit of variables. + + Parameters + ---------- + var_name + name of variable + dataset + name of dataset + + Returns + ------- + str + unit of variable + """ + if var_name == "psl": + return "[hPa]" + if var_name == "pr" and dataset == "ERA5": + return "[m]" + if var_name == "pr": + return "[kg m-2 s-1]" + if var_name == "tas": + return "[K]" + return "" + + +def plot_maps( + wt: np.array, + cfg: dict, + cube: iris.cube.Cube, + data_info: dict, + mode: str, +): + """Plot maps of means, std and anomalies. + + Parameters + ---------- + wt + WT array + cfg + Configuration dictionary from recipe + cube + Data to be plotted + data_info + Dictionary with info to dataset + mode + Statistics that is used + """ + var_name = data_info.get("var") + + logger.info( + "Plotting %s %s %s for %s %s", + data_info.get("dataset"), + var_name, + mode, + data_info.get("wt_string"), + wt, + ) + + ax = plt.axes(projection=ccrs.PlateCarree()) + + cmap = get_colormap(var_name if var_name != "pr" else "prcp") + + title = prepare_plot_title(data_info, wt, var_name, mode) + plt.title(title) + unit = get_unit(var_name, data_info.get("dataset")) + + im = iplt.contourf(cube if var_name != "psl" else cube / 100, cmap=cmap) + cb = plt.colorbar(im) + cb.ax.tick_params(labelsize=8) + cb.set_label(label=f"{var_name} {mode} {unit}") + + set_gridlines(ax) + + ax.set_extent([-15, 20, 27.5, 62.5]) + + ax.coastlines() + ax.add_feature(cfeature.BORDERS, linestyle=":") + + if not Path(f"{cfg.get('plot_dir')}/{mode}").exists(): + Path(f"{cfg.get('plot_dir')}/{mode}").mkdir( + parents=True, + exist_ok=True, + ) + + plt.savefig( + f"{cfg.get('plot_dir')}/{mode}/{data_info.get('wt_string')}_{wt}{get_driver(data_info)}_{data_info.get('dataset')}_{data_info.get('ensemble', '')}" + f"_{var_name}_{mode}_{data_info.get('timerange')}.png", + ) + plt.savefig( + f"{cfg.get('plot_dir')}/{mode}/{data_info.get('wt_string')}_{wt}{get_driver(data_info)}_{data_info.get('dataset')}_{data_info.get('ensemble', '')}_" + f"{var_name}_{mode}_{data_info.get('timerange')}.pdf", + ) + plt.close() + + # log provenance + ancestors = [ + f"{data_info.get('preproc_path')}", + f"{cfg.get('work_dir')}/ERA5.nc", + ] + provenance_record = get_provenance_record( + f"{var_name} {mode} for {data_info.get('wt_string')}", + ancestors, + [var_name], + ["map"], + [mode], + ) + + local_path = f"{cfg.get('plot_dir')}/{mode}" + + log_provenance( + f"{local_path}/{data_info.get('wt_string')}_{wt}{get_driver(data_info)}_" + f"{data_info.get('dataset')}_{data_info.get('ensemble', '')}" + f"_{var_name}_{mode}_{data_info.get('timerange')}.png", + cfg, + provenance_record, + ) + + +def plot_corr_rmse_heatmaps( + cfg: dict, + pattern_correlation_matrix: np.array, + rmse_matrix: np.array, + dataset: str, + timerange: str, +): + """Plot correlation and rmse heatmaps. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + pattern_correlation_matrix + Pattern correlation matrix + rmse_matrix + RMSE matrix + dataset + Name of dataset + timerange + Time range for the calculation + """ + output_path = f"{cfg.get('plot_dir')}/heatmaps" + + if not Path(f"{output_path}").exists(): + Path(f"{output_path}").mkdir(parents=True, exist_ok=True) + + labels = np.arange(1, 28) + + mask = np.zeros_like(pattern_correlation_matrix) + mask[np.triu_indices_from(mask)] = True + with sns.axes_style("white"): + plt.figure(figsize=(10, 10)) + plt.title(f"Correlation Matrix , {dataset}, {timerange}") + levels = np.linspace( + np.min(pattern_correlation_matrix), + np.max(pattern_correlation_matrix), + 9, + ) + ax = sns.heatmap( + pattern_correlation_matrix, + mask=mask, + square=True, + annot=True, + annot_kws={"size": 6}, + cmap="seismic", + xticklabels=labels, + yticklabels=labels, + cbar_kws={"ticks": levels, "shrink": 0.8, "format": "%.2f"}, + ) + ax.set_xlabel("lwt", fontsize=8) + ax.set_xticklabels(ax.get_xticklabels(), rotation=0) + ax.set_yticklabels(ax.get_yticklabels(), rotation=0) + ax.set_ylabel("lwt", fontsize=8) + plt.tight_layout() + plt.savefig( + f"{output_path}/correlation_matrix_{dataset}_{timerange}.png", + ) + plt.savefig( + f"{output_path}/correlation_matrix_{dataset}_{timerange}.pdf", + ) + plt.close() + + mask = np.zeros_like(rmse_matrix) + mask[np.triu_indices_from(mask)] = True + with sns.axes_style("white"): + plt.figure(figsize=(10, 10)) + plt.title(f"RMSE Matrix, {dataset}, {timerange}") + levels = np.linspace(np.min(rmse_matrix), np.max(rmse_matrix), 9) + ax = sns.heatmap( + rmse_matrix, + mask=mask, + square=True, + annot=True, + annot_kws={"size": 5}, + cmap="seismic", + xticklabels=labels, + yticklabels=labels, + cbar_kws={"ticks": levels, "shrink": 0.8, "format": "%.2f"}, + ) + ax.set_xlabel("lwt", fontsize=8) + ax.set_xticklabels(ax.get_xticklabels(), rotation=0) + ax.set_yticklabels(ax.get_yticklabels(), rotation=0) + ax.set_ylabel("lwt", fontsize=8) + plt.tight_layout() + plt.savefig(f"{output_path}/rmse_matrix_{dataset}_{timerange}.png") + plt.savefig(f"{output_path}/rmse_matrix_{dataset}_{timerange}.pdf") + plt.close() + # log provenance + # rmse matrix + ancestors = [ + f"{cfg.get('work_dir')}/ERA5.nc", + ] + provenance_record = get_provenance_record( + "rmse matrix", + ancestors, + ["rmse"], + ["other"], + ) + + local_path = f"{cfg.get('plot_dir')}/" + + log_provenance( + f"{local_path}/rmse_matrix_{dataset}_{timerange}.png", + cfg, + provenance_record, + ) + # correlation matrix + provenance_record = get_provenance_record( + "correlation matrix", + ancestors, + ["correlation"], + ["other"], + ) + + log_provenance( + f"{local_path}/correlation_matrix_{dataset}_{timerange}.png", + cfg, + provenance_record, + ) + + +def get_colormap(colormap_string: str) -> ListedColormap: + """Get colormaps for plottings. + + Parameters + ---------- + colormap_string + string to identify colormap + + Returns + ------- + ListedColormap + Colormap for the specified variable + """ + misc_seq_2_disc = [ + (230 / 255, 240 / 255, 240 / 255), + (182 / 255, 217 / 255, 228 / 255), + (142 / 255, 192 / 255, 226 / 255), + (118 / 255, 163 / 255, 228 / 255), + (116 / 255, 130 / 255, 222 / 255), + (121 / 255, 97 / 255, 199 / 255), + (118 / 255, 66 / 255, 164 / 255), + (107 / 255, 40 / 255, 121 / 255), + (86 / 255, 22 / 255, 75 / 255), + (54 / 255, 14 / 255, 36 / 255), + ] + + temp_seq_disc = [ + (254 / 255, 254 / 255, 203 / 255), + (251 / 255, 235 / 255, 153 / 255), + (244 / 255, 204 / 255, 104 / 255), + (235 / 255, 167 / 255, 84 / 255), + (228 / 255, 134 / 255, 80 / 255), + (209 / 255, 98 / 255, 76 / 255), + (164 / 255, 70 / 255, 66 / 255), + (114 / 255, 55 / 255, 46 / 255), + (66 / 255, 40 / 255, 24 / 255), + (25 / 255, 25 / 255, 0 / 255), + ] + + prec_seq_disc = [ + (255 / 255, 255 / 255, 229 / 255), + (217 / 255, 235 / 255, 213 / 255), + (180 / 255, 216 / 255, 197 / 255), + (142 / 255, 197 / 255, 181 / 255), + (105 / 255, 177 / 255, 165 / 255), + (67 / 255, 158 / 255, 149 / 255), + (44 / 255, 135 / 255, 127 / 255), + (29 / 255, 110 / 255, 100 / 255), + (14 / 255, 85 / 255, 74 / 255), + (0 / 255, 60 / 255, 48 / 255), + ] + + if colormap_string == "psl": + return ListedColormap(misc_seq_2_disc) + if colormap_string == "prcp": + return ListedColormap(prec_seq_disc) + if colormap_string == "temp": + return ListedColormap(temp_seq_disc) + + return None diff --git a/esmvaltool/diag_scripts/weathertyping/weathertyping.py b/esmvaltool/diag_scripts/weathertyping/weathertyping.py new file mode 100644 index 0000000000..0f8da4d277 --- /dev/null +++ b/esmvaltool/diag_scripts/weathertyping/weathertyping.py @@ -0,0 +1,409 @@ +"""Script for calculating Lamb weathertypes. + +It also plots the means and seasonal occurrence of the +weathertypes, and offers the option to calculate simplified +weathertypes based on precipitation patterns. +""" + +import iris + +from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic +from esmvaltool.diag_scripts.weathertyping.calc_utils import ( + calc_lwt_model, + calc_lwt_slwt_model, + calc_slwt_obs, + plot_means, + wt_algorithm, +) +from esmvaltool.diag_scripts.weathertyping.plot_utils import ( + plot_seasonal_occurrence, +) +from esmvaltool.diag_scripts.weathertyping.wt_utils import ( + combine_wt_to_file, + get_ancestors_era5_eobs, + get_looping_dict, + get_model_output_filepath, + get_preproc_lists_ensemble, + get_provenance_record, + load_wt_files, + load_wt_preprocessors, + log_provenance, + run_predefined_slwt, + write_lwt_to_file, +) + + +def process_models_automatic_slwt( + cfg: dict, + dataset_vars: list, + data_info: dict, +): + """Process model data for calculating Lamb and simplified weathertypes. + + Parameters + ---------- + cfg + Nested dictionary of metadata. + dataset_vars + List of variable dictionaries for a specific dataset. + data_info + Dictionary holding dataset information. + """ + for ensemble_var in dataset_vars: + if ensemble_var.get("preprocessor") == "weathertype_preproc": + data_info["ensemble"] = ensemble_var.get("ensemble", "") + data_info["driver"] = ensemble_var.get("driver", "") + + if data_info["driver"] != "": + data_info["driver"] = "_" + data_info["driver"] + + wt_preproc = iris.load_cube(ensemble_var.get("filename")) + + output_file_path, preproc_path = get_model_output_filepath( + data_info["dataset"], + ensemble_var, + ) + + data_info["output_file_path"] = output_file_path + data_info["preproc_path"] = preproc_path + + # calculate weathertypes + calc_lwt_slwt_model( + cfg, + wt_preproc, + data_info, + cfg.get("predefined_slwt"), + ) + # plot means + if cfg.get("plotting", False): + # load wt files + wt_cube_path = ( + f"{cfg.get('work_dir')}/{output_file_path}" + f"/{data_info['dataset']}" + f"{data_info['driver']}_" + f"{data_info['timerange']}.nc" + ) + + wt_cubes = load_wt_files(wt_cube_path) + + var_dict = { + f"{ensemble_var.get('short_name')}": get_preproc_lists_ensemble( + ensemble_var, + ), + } + for var_name, var_data in var_dict.items(): + data_info["var"] = var_name + data_info["preproc_path"] = var_data[1] + + plot_means(cfg, var_data[0], wt_cubes, data_info) + plot_seasonal_occurrence( + cfg, + wt_cubes, + data_info, + wt_cube_path, + ) + + +def process_era5_automatic_slwt( + data_info: dict, + preproc_variables_dict: dict, + cfg: dict, + dataset_vars: list, +): + """Process ERA5 data for calculating Lamb and simplified weathertypes. + + Parameters + ---------- + data_info + Dictionary holding dataset information. + preproc_variables_dict + Dictionary holding preprocessed variables for all datasets. + cfg + Nested dictionary of metadata. + dataset_vars + List of variable dictionaries for a specific dataset. + """ + wt_preproc, wt_preproc_prcp, wt_preproc_prcp_eobs = load_wt_preprocessors( + data_info["dataset"], + preproc_variables_dict, + ) + + # calculate lwt + lwt = wt_algorithm(wt_preproc, data_info["dataset"]) + + era5_ancestors, eobs_ancestors = get_ancestors_era5_eobs( + data_info["dataset"], + preproc_variables_dict, + ) + + # calculate simplified lwt based on precipitation + # patterns or use predefined_slwt + if not cfg.get("predefined_slwt"): + slwt_era5 = calc_slwt_obs( + cfg, + lwt, + wt_preproc_prcp, + data_info["dataset"], + era5_ancestors, + data_info["timerange"], + ) + slwt_eobs = calc_slwt_obs( + cfg, + lwt, + wt_preproc_prcp_eobs, + "E-OBS", + eobs_ancestors, + data_info["timerange"], + ) + else: + slwt_era5, slwt_eobs = run_predefined_slwt( + cfg.get("work_dir"), + data_info["dataset"], + lwt, + cfg.get("predefined_slwt"), + ) + + # write to file + wt_list = [lwt, slwt_era5, slwt_eobs] + combine_wt_to_file(cfg, wt_list, wt_preproc, data_info["dataset"]) + + # load weathertype files as cubes + wt_cube_path = f"{cfg.get('work_dir')}/{data_info['dataset']}.nc" + wt_cubes = load_wt_files(wt_cube_path) + + if cfg.get("plotting", False): + var_dict = get_looping_dict( + dataset_vars, + ) # dataset_vars is list of variables for dataset dataset_name + # plot means + for var_name, var_data in var_dict.items(): + data_info["var"] = var_name + data_info["preproc_path"] = var_data[1] + + plot_means(cfg, var_data[0], wt_cubes, data_info) + plot_seasonal_occurrence(cfg, wt_cubes, data_info, wt_cube_path) + + +def run_automatic_slwt(cfg: dict): + """Run the automated calculation for simplified weathertypes and write to file, and plot the means and seasonal occurrence of the weathertypes. + + Parameters + ---------- + cfg + Nested dictionary of metadata + """ + preproc_variables_dict = group_metadata( + cfg.get("input_data").values(), + "dataset", + ) + for dataset_name, dataset_vars in preproc_variables_dict.items(): + data_info = { + "timerange": dataset_vars[0].get("timerange").replace("/", "-"), + "dataset": dataset_name, + } + if dataset_name == "ERA5": + process_era5_automatic_slwt( + data_info, + preproc_variables_dict, + cfg, + dataset_vars, + ) + else: + if data_info["dataset"] == "E-OBS": + continue + process_models_automatic_slwt(cfg, dataset_vars, data_info) + + +def process_era5_lwt( + preproc_variables_dict: dict, + cfg: dict, + dataset_vars: list, + data_info: dict, +): + """Process ERA5 data for calculating Lamb weathertypes. + + Parameters + ---------- + preproc_variables_dict + Dictionary holding preprocessed variables for all datasets. + cfg + Nested dictionary of metadata. + dataset_vars + List of variable dictionaries for a specific dataset. + data_info + Dictionary holding dataset information. + """ + wt_preproc, _, _ = load_wt_preprocessors( + data_info["dataset"], + preproc_variables_dict, + ) + + # calculate lwt + lwt = wt_algorithm(wt_preproc, data_info["dataset"]) + + era5_ancestors, eobs_ancestors = get_ancestors_era5_eobs( + data_info["dataset"], + preproc_variables_dict, + ) + + ancestors = [era5_ancestors, eobs_ancestors] + + provenance_record = get_provenance_record( + "Lamb Weathertypes", + ancestors, + ["Lamb Weathertypes"], + ) + + log_provenance(f"{cfg.get('work_dir')}/lwt_era5", cfg, provenance_record) + + # write only lwt to file + write_lwt_to_file(cfg, lwt, wt_preproc, data_info["dataset"]) + + # load wt files + wt_cube_path = f"{cfg.get('work_dir')}/{data_info['dataset']}.nc" + wt_cubes = load_wt_files( + wt_cube_path, + mode="lwt", + ) + + if cfg.get("plotting", False): + var_dict = get_looping_dict( + dataset_vars, + ) # dataset_vars is list of variables for dataset dataset_name + # plot means + for var_name, var_data in var_dict.items(): + data_info["var"] = var_name + data_info["preproc_path"] = var_data[1] + + plot_means(cfg, var_data[0], wt_cubes, data_info, "lwt") + plot_seasonal_occurrence(cfg, wt_cubes, data_info, wt_cube_path) + + +def process_models_lwt(cfg: dict, dataset_vars: list, data_info: dict): + """Process model data for calculating Lamb weathertypes. + + Parameters + ---------- + cfg + Nested dictionary of metadata. + dataset_vars + List of variable dictionaries for a specific dataset. + data_info + Dictionary holding dataset information. + """ + for ensemble_var in dataset_vars: + if ensemble_var.get("preprocessor") == "weathertype_preproc": + data_info["ensemble"] = ensemble_var.get("ensemble", "") + data_info["driver"] = ensemble_var.get("driver", "") + + if data_info["driver"] != "": + data_info["driver"] = "_" + {data_info["driver"]} + + wt_preproc = iris.load_cube(dataset_vars[0].get("filename")) + + output_file_path, preproc_path = get_model_output_filepath( + data_info["dataset"], + dataset_vars, + ) + + data_info["output_file_path"] = output_file_path + data_info["preproc_path"] = preproc_path + + # calculate weathertypes + calc_lwt_model(cfg, wt_preproc, data_info) + + # load wt files + wt_cube_path = ( + f"{cfg.get('work_dir')}/{output_file_path}" + f"/{data_info['dataset']}" + f"{data_info['driver']}" + f"_{data_info['ensemble']}_" + f"{data_info['timerange']}.nc" + ) + + wt_cubes = load_wt_files( + wt_cube_path, + mode="lwt", + ) + + var_dict = get_looping_dict( + dataset_vars, + ) # dataset_vars is list of variables for dataset_name + + if cfg.get("plotting", False): + # plot means + for var_name, var_data in var_dict.items(): + data_info["var"] = var_name + data_info["preproc_path"] = var_data[1] + + plot_means( + cfg, + var_data[0], + wt_cubes, + data_info, + "lwt", + ) + plot_seasonal_occurrence( + cfg, + wt_cubes, + data_info, + wt_cube_path, + ) + + +def run_lwt(cfg: dict): + """Run calculation of weathertypes. Write to file, and plot the means of psl, tas, and pr for each weathertype. Plot seasonal occurrence of the weathertypes. + + Parameters + ---------- + cfg + Nested dictionary of metadata. + """ + preproc_variables_dict = group_metadata( + cfg.get("input_data").values(), + "dataset", + ) + for dataset_name, dataset_vars in preproc_variables_dict.items(): + data_info = { + "timerange": dataset_vars[0].get("timerange").replace("/", "-"), + "dataset": dataset_name, + } + + if dataset_name == "ERA5": + process_era5_lwt( + preproc_variables_dict, + cfg, + dataset_vars, + data_info, + ) + else: + if data_info["dataset"] == "E-OBS": + continue + process_models_lwt(cfg, dataset_vars, data_info) + + +def run_weathertyping_diagnostic(cfg: dict): + """Run the weathertyping diagnostic. + + Parameters + ---------- + cfg + Nested dictionary of metadata + """ + automatic_slwt = cfg.get("automatic_slwt") + + # check if user wants to calculate simplified weathertypes automatically + # for that, automatic_slwt must be true + # additionally, if a predefined_slwt is given, those will be used + if automatic_slwt: + run_automatic_slwt(cfg) + # if automatic_slwt is false, and predefined_slwt is false, + # just write selected pairs to file + else: + run_lwt(cfg) + + +if __name__ == "__main__": + with run_diagnostic() as config: + # main function for running the diagnostic + run_weathertyping_diagnostic(config) diff --git a/esmvaltool/diag_scripts/weathertyping/wt_utils.py b/esmvaltool/diag_scripts/weathertyping/wt_utils.py new file mode 100644 index 0000000000..4f5a8e99ca --- /dev/null +++ b/esmvaltool/diag_scripts/weathertyping/wt_utils.py @@ -0,0 +1,632 @@ +"""Utility functions for weathertyping script.""" + +import json +import logging +import warnings +from pathlib import Path + +import iris +import iris.analysis.cartography +import numpy as np +import pandas as pd + +from esmvaltool.diag_scripts.shared import ProvenanceLogger + +iris.FUTURE.datum_support = True + +logger = logging.getLogger(Path(__file__).name) + +# Ignoring a warning that is produced when selecting timesteps of a weathertype +warnings.filterwarnings("ignore", ".*Collapsing a non-contiguous coordinate*") + + +def get_driver(data_info: dict) -> str: + """Get driving model name and string for further use. + + Parameters + ---------- + data_info + Data information dictionary. + + Returns + ------- + str: + Driver string with leading underscore or empty string. + """ + return data_info.get("driver", "") + + +def load_wt_preprocessors(dataset: str, preproc_variables_dict: dict) -> tuple: + """Load preprocessor cubes for calculating Lamb weathertypes. + + Parameters + ---------- + dataset + Name of dataset + preproc_variables_dict + Dictionary with info on preprocessor variables + + Returns + ------- + tuple + Preprocessor cubes for weathertyping + """ + wt_preproc = iris.load_cube( + preproc_variables_dict.get(dataset)[0].get("filename"), + ) + # try and except block to catch missing precipitation preprocessor + # if it is not supplied, predefined_slwt has to be given by user + try: + wt_preproc_prcp = iris.load_cube( + preproc_variables_dict.get(dataset)[1].get("filename"), + ) + wt_preproc_prcp_eobs = iris.load_cube( + preproc_variables_dict.get("E-OBS")[0].get("filename"), + ) + except (IndexError, KeyError, TypeError): + logger.info( + "ERA5 precipitation preprocessor not found for automatic slwt.", + ) + wt_preproc_prcp = None + wt_preproc_prcp_eobs = None + + return wt_preproc, wt_preproc_prcp, wt_preproc_prcp_eobs + + +def get_ancestors_era5_eobs( + dataset: str, + preproc_variables_dict: dict, +) -> tuple: + """Get ancestors for observational data. + + Parameters + ---------- + dataset + Name of dataset + preproc_variables_dict + Dictionary with info on preprocessor variables + + Returns + ------- + tuple + Lists of ERA5 and E-OBS ancestors + """ + era5_ancestors = [ + preproc_variables_dict.get(dataset)[0].get("filename"), + ] + try: + era5_ancestors.append( + preproc_variables_dict.get(dataset)[1].get("filename"), + ) + eobs_ancestors = [ + preproc_variables_dict.get(dataset)[0].get("filename"), + preproc_variables_dict.get("E-OBS")[0].get("filename"), + ] + except (IndexError, KeyError, TypeError): + logger.info( + "No ancestors for ERA5 and E-OBS precipitation preprocessor.", + ) + eobs_ancestors = [] + + return era5_ancestors, eobs_ancestors + + +def get_model_output_filepath(dataset: str, data_info: list) -> tuple: + """Get output filepaths for models. + + Parameters + ---------- + dataset + Name of dataset + data_info + Model variables + + Returns + ------- + tuple + Output filepath and preprocessor path for future referencing. + """ + timerange = data_info.get("timerange").replace("/", "-") + experiment = data_info.get("exp") + ensemble = data_info.get("ensemble") + driver = data_info.get("driver", "") + + out_path = f"{dataset}/{driver}/{experiment}/{ensemble}/{timerange}" + preproc_path = data_info.get("filename") + + return out_path, preproc_path + + +def get_preproc_lists(preproc_vars: list): + """Put preprocessors and paths into lists for further use. + + Parameters + ---------- + preproc_vars + List of preprocessor variables. + + Returns + ------- + tuple + Preprocessor cubes and paths. + """ + preproc_path_psl = preproc_vars[-3].get("filename") + preproc_path_prcp = preproc_vars[-2].get("filename") + preproc_path_tas = preproc_vars[-1].get("filename") + + mean_preproc_psl = iris.load_cube(preproc_vars[-3].get("filename")) + mean_preproc_prcp = iris.load_cube(preproc_vars[-2].get("filename")) + mean_preproc_tas = iris.load_cube(preproc_vars[-1].get("filename")) + + preproc_list = [mean_preproc_psl, mean_preproc_prcp, mean_preproc_tas] + preproc_path_list = [preproc_path_psl, preproc_path_prcp, preproc_path_tas] + + return preproc_list, preproc_path_list + + +def get_preproc_lists_ensemble(preproc_vars: list): + """Put preprocessors and paths into lists for further use. + + Parameters + ---------- + preproc_vars + List of preprocessor variables. + + Returns + ------- + tuple + Preprocessor cubes and paths. + """ + preproc_path = preproc_vars.get("filename") + + preproc = iris.load_cube(preproc_vars.get("filename")) + + return preproc, preproc_path + + +def get_looping_dict(preproc_vars: list): + """Put variable preprocessors into dict for looping. + + Parameters + ---------- + preproc_vars + List of preprocessor variables. + + Returns + ------- + dict + Dictionary of preprocessor cubes and paths. + """ + preproc, preproc_path = get_preproc_lists(preproc_vars) + + return { + "psl": [preproc[0], preproc_path[0]], + "pr": [preproc[1], preproc_path[1]], + "tas": [preproc[2], preproc_path[2]], + } + + +def load_wt_files(path: str, mode="slwt"): + """Load wt files. + + Parameters + ---------- + path + Path of wt file + mode + Type of weathertype to load. Defaults to "slwt". + + Returns + ------- + list + List of weathertype cubes. + """ + if mode == "slwt": + lwt_cube = iris.load_cube(path, "lwt") + slwt_era5_cube = iris.load_cube(path, "slwt_era5") + slwt_eobs_cube = iris.load_cube(path, "slwt_eobs") + wt_cubes = [lwt_cube, slwt_era5_cube, slwt_eobs_cube] + elif mode == "lwt": + lwt_cube = iris.load_cube(path, "lwt") + wt_cubes = [lwt_cube] + else: + e = "Mode not recognized. Use 'slwt' or 'lwt'." + raise ValueError(e) + + return wt_cubes + + +def get_provenance_record( + caption: str, + ancestors: list, + long_names: list, + plot_types: list | None = None, + statistics: list | None = None, +) -> dict: + """Get provenance record. + + Parameters + ---------- + caption + Caption of plot + ancestors + List of ancestor plots + long_names + List of variable long names + plot_types + Type of plot + statistics + Types of statistics used + + Returns + ------- + dict + Provenance record + """ + record = { + "caption": caption, + "domains": ["reg"], + "authors": ["jury_martin", "kroissenbrunner_thomas"], + "references": ["maraun21jgr", "jones93ijc"], + "projects": ["preval"], + "long_names": long_names, + "ancestors": ancestors, + } + if plot_types: + record["plot_types"] = plot_types + if statistics: + record["statistics"] = statistics + return record + + +def log_provenance(filename: str, cfg: dict, provenance_record: dict): + """Log provenance. + + Parameters + ---------- + filename: + Filename of xml file + cfg + Configuration dictionary + provenance_record + Provenance record dictionary + """ + with ProvenanceLogger(cfg) as provenance_logger: + provenance_logger.log(filename, provenance_record) + + logger.info( + "Provenance added to %s", + f"{cfg['run_dir']}/diagnostic_provenance.yml", + ) + + +def turn_list_to_mapping_dict(list_: list) -> dict: + """Turn list of combined WT to a dictionary for further use. + + Parameters + ---------- + list_ + List of combined WTs + + Returns + ------- + dict + Dictionary of combined WTs + """ + result_dict = {} + + for i, s in enumerate(list_): + for elem in s: + if elem not in result_dict: + result_dict[elem] = i + 1 + + return result_dict + + +def get_mapping_dict(selected_pairs: list) -> dict: + """Get mapping dictionary from list of selected pairs. + + Parameters + ---------- + selected_pairs + List of selected weathertype pairs + + Returns + ------- + dict + Mapping dictionary + """ + # selected pairs is of form [(wt1, wt2), (wt3, wt4), ...] + sets_list = [set(i) for i in selected_pairs if i] + + # 3 loops to check all sets against each other until no more merges happen + # merged_flag indicates if a merge happened in the last full pass -> we need to restart then + merged_flag = True + + while merged_flag: + merged_flag = False + i = 0 + + # check each set against each other + while i < len(sets_list): + set1 = sets_list[i] + j = i + 1 + + while j < len(sets_list): + set2 = sets_list[j] + + if set1.isdisjoint(set2) is False: + # this means there is an overlap -> merge sets + set1.update(set2) + sets_list.pop(j) + merged_flag = True + else: + j += 1 + + i += 1 + + return turn_list_to_mapping_dict(sets_list) + + +def write_mapping_dict(work_dir: str, dataset: str, mapping_dict: dict): + """Write mapping dictionary to file. + + Parameters + ---------- + work_dir + Current working directory + dataset + Dataset name + mapping_dict + Mapping dictionary in {lwt: slwt, ...} format + """ + mapping_dict_reformat = convert_dict(mapping_dict) + + with Path( + f"{work_dir}/wt_mapping_dict_{dataset}.json", + ).open("w", encoding="utf-8") as file: + json.dump(mapping_dict_reformat, file) + + +def convert_dict(dict_: dict) -> dict: + """Convert dictionary from {lwt: slwt, ...} format to {slwt: [lwt1, lwt2], ...}. + + Parameters + ---------- + dict_ + Mapping dictionary to be converted + + Returns + ------- + dict + Converted dictionary + """ + new_dict = {} + for dataset, value in dict_.items(): + if value not in new_dict: + new_dict[value] = [] + new_dict[value].append(dataset) + return new_dict + + +def reverse_convert_dict(originial_dict: dict) -> dict: + """Convert mapping dictionary. + + From {slwt: [lwt1, lwt2], ...} format to {lwt: slwt, ...}. + + Parameters + ---------- + original_dict + Dict in the {slwt: [lwt1, lwt2], ...} format + + Returns + ------- + dict + Dict in the format {lwt: slwt, ...} + """ + new_dict = {} + for key, value_list in originial_dict.items(): + for original_key in value_list: + new_dict[original_key] = key + return new_dict + + +def write_corr_rmse_to_csv( + cfg: dict, + pattern_correlation_matrix: np.array, + rmse_matrix: np.array, + dataset: str, +): + """Write correlation and rmse matrix to csv files. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + pattern_correlation_matrix + Correlation matrix + rmse_matrix + RMSE matrix + dataset + Name of dataset + """ + logger.info("Writing corr and rsme matrices for %s", dataset) + + work_dir = cfg.get("work_dir") + + df_corr = pd.DataFrame(pattern_correlation_matrix) + df_corr.index = range(1, len(df_corr) + 1) + df_corr.columns = range(1, len(df_corr.columns) + 1) + df_corr.to_csv( + f"{work_dir}/correlation_matrix_{dataset}.csv", + index_label="Index", + ) + + df_rmse = pd.DataFrame(rmse_matrix) + df_rmse.index = range(1, len(df_rmse) + 1) + df_rmse.columns = range(1, len(df_rmse.columns) + 1) + df_rmse.to_csv( + f"{work_dir}/rmse_matrix_{dataset}.csv", + index_label="Index", + ) + + +def run_predefined_slwt( + work_dir: str, + dataset_name: str, + lwt: np.array, + predefined_slwt: dict, +): + """Run predefined slwt mapping. + + Parameters + ---------- + work_dir + Working directory to save mapping dict + dataset_name + Name of dataset + lwt + lwt array + predefined_slwt + Mapping dictionary in {lwt: slwt, ...} format + + Returns + ------- + np.array + slwt_era5 array + np.array + slwt_eobs array + """ + predefined_slwt = check_mapping_dict_format(predefined_slwt) + write_mapping_dict(work_dir, dataset_name, predefined_slwt) + write_mapping_dict(work_dir, "E-OBS", predefined_slwt) + slwt_era5 = map_lwt_to_slwt(lwt, predefined_slwt) + slwt_eobs = map_lwt_to_slwt(lwt, predefined_slwt) + + return slwt_era5, slwt_eobs + + +def combine_wt_to_file( + cfg: dict, + wt_list: list, + cube: iris.cube.Cube, + file_name: str, +): + """Combine lwt and slwt arrays to one file. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + wt_list + List of weathertype arrays + cube + Cube of data to keep time coordinate + file_name + Name of output file + """ + lwt = wt_list[0] + slwt_era5 = wt_list[1] + slwt_eobs = wt_list[2] + + logger.info("Writing weathertypes to %s", file_name) + + tcoord = cube.coord("time") + time_points = tcoord.units.num2date(tcoord.points) + + write_path = cfg.get("work_dir") + + wt_cube = iris.cube.CubeList() + wt_cube.append(iris.cube.Cube(lwt, long_name="lwt")) + wt_cube.append(iris.cube.Cube(slwt_era5, long_name="slwt_era5")) + wt_cube.append(iris.cube.Cube(slwt_eobs, long_name="slwt_eobs")) + + wt_cube[0].add_dim_coord(tcoord, 0) + wt_cube[1].add_dim_coord(tcoord, 0) + wt_cube[2].add_dim_coord(tcoord, 0) + + iris.save(wt_cube, f"{write_path}/{file_name}.nc") + + # write to csv file + d = { + "date": time_points[:], + "lwt": np.int8(lwt), + "slwt_era5": np.int8(slwt_era5), + "slwt_eobs": np.int8(slwt_eobs), + } + df = pd.DataFrame(data=d) + df.to_csv(write_path + f"/{file_name}.csv", index=False) + + +def write_lwt_to_file( + cfg: dict, + lwt: np.array, + cube: iris.cube.Cube, + file_name: str, +): + """Write only lwt to file. + + Parameters + ---------- + cfg + Configuration dictionary from recipe + lwt + lwt array + cube + Cube of data to keep time coordinate + file_name + Name of output file + """ + logger.info("Writing Lamb Weathertype to %s", file_name) + + tcoord = cube.coord("time") + time_points = tcoord.units.num2date(tcoord.points) + + write_path = cfg.get("work_dir") + + wt_cube = iris.cube.CubeList() + wt_cube.append(iris.cube.Cube(np.int8(lwt), long_name="lwt")) + + wt_cube[0].add_dim_coord(tcoord, 0) + iris.save(wt_cube, f"{write_path}/{file_name}.nc") + + # write to csv file + d = {"date": time_points[:], "lwt": np.int8(lwt)} + df = pd.DataFrame(data=d) + df.to_csv(write_path + f"/{file_name}.csv", index=False) + + +def map_lwt_to_slwt(lwt: np.array, mapping_dict: dict) -> np.array: + """Map lwt array to slwt array. + + Parameters + ---------- + lwt + lwt array + mapping_dict + Mapping dictionary in {lwt: slwt, ...} format + + Returns + ------- + np.array + array of slwt + """ + return np.array([np.int8(mapping_dict.get(value, 0)) for value in lwt]) + + +def check_mapping_dict_format(mapping_dict: dict) -> dict: + """Check format of mapping dict and return in {lwt: slwt, ...} format. + + Parameters + ---------- + mapping_dict + mapping dict in any format + + Returns + ------- + dict + mapping dict in {lwt: slwt, ...} format + """ + if isinstance(mapping_dict[next(iter(mapping_dict))], list): + return reverse_convert_dict(mapping_dict) + return mapping_dict diff --git a/esmvaltool/recipes/recipe_weathertyping_CMIP6.yml b/esmvaltool/recipes/recipe_weathertyping_CMIP6.yml new file mode 100644 index 0000000000..c745ff5b3e --- /dev/null +++ b/esmvaltool/recipes/recipe_weathertyping_CMIP6.yml @@ -0,0 +1,305 @@ +documentation: + title: Weathertyping algorithm + + description: | + A diagnostic to calculate Lamb weathertypes over a given region. Furthermore, + correlations between weathertypes and precipitation patterns over a given area can be calculated + and 'combined' or 'simplified' weathertypes can be derived. Additionally, mean fields, as well as + anomalies and standard deviations can be plotted. + + authors: + - jury_martin + - kroissenbrunner_thomas + + maintainer: + - kroissenbrunner_thomas + + projects: + - preval + + references: + - maraun21jgr + - jones93ijc + +datasets_psl_CMIP6: &datasets_psl_cmip6 + - {institute: AS-RCEC, dataset: TaiESM1, ensemble: r1i1p1f1, grid: gn, esgf_version: v20200626} + - {institute: AWI, dataset: AWI-ESM-1-1-LR, ensemble: r1i1p1f1, grid: gn, esgf_version: v20200212} + - {institute: BCC, dataset: BCC-CSM2-MR, ensemble: r1i1p1f1, grid: gn, esgf_version: v20181126} + - {institute: BCC, dataset: BCC-ESM1, ensemble: r1i1p1f1, grid: gn, esgf_version: v20181220} + - {institute: CAS, dataset: FGOALS-f3-L, ensemble: r1i1p1f1, grid: gr, esgf_version: v20191019} + - {institute: CAS, dataset: FGOALS-g3, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190826} + - {institute: CCCma, dataset: CanESM5, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190429} + - {institute: CMCC, dataset: CMCC-CM2-HR4, ensemble: r1i1p1f1, grid: gn, esgf_version: v20200904} + - {institute: CMCC, dataset: CMCC-CM2-SR5, ensemble: r1i1p1f1, grid: gn, esgf_version: v20200616} + - {institute: CMCC, dataset: CMCC-ESM2, ensemble: r1i1p1f1, grid: gn, esgf_version: v20210114} + - {institute: CNRM-CERFACS, dataset: CNRM-CM6-1, ensemble: r1i1p1f2, grid: gr, esgf_version: v20180917} + - {institute: CSIRO-ARCCSS, dataset: ACCESS-CM2, ensemble: r1i1p1f1, grid: gn, esgf_version: v20191108} + - {institute: CSIRO, dataset: ACCESS-ESM1-5, ensemble: r1i1p1f1, grid: gn, esgf_version: v20191115} + - {institute: EC-Earth-Consortium, dataset: EC-Earth3-AerChem, ensemble: r1i1p1f1, grid: gr, esgf_version: v20200624} + - {institute: EC-Earth-Consortium, dataset: EC-Earth3-CC, ensemble: r1i1p1f1, grid: gr, esgf_version: v20210113} + - {institute: EC-Earth-Consortium, dataset: EC-Earth3-Veg-LR, ensemble: r1i1p1f1, grid: gr, esgf_version: v20200217} + - {institute: EC-Earth-Consortium, dataset: EC-Earth3-Veg, ensemble: r1i1p1f1, grid: gr, esgf_version: v20211207} + - {institute: EC-Earth-Consortium, dataset: EC-Earth3, ensemble: r1i1p1f1, grid: gr, esgf_version: v20200310} + - {institute: HAMMOZ-Consortium, dataset: MPI-ESM-1-2-HAM, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190627} + - {institute: INM, dataset: INM-CM4-8, ensemble: r1i1p1f1, grid: gr1, esgf_version: v20190530} + - {institute: INM, dataset: INM-CM5-0, ensemble: r1i1p1f1, grid: gr1, esgf_version: v20190610} + - {institute: IPSL, dataset: IPSL-CM5A2-INCA, ensemble: r1i1p1f1, grid: gr, esgf_version: v20200729} + - {institute: IPSL, dataset: IPSL-CM6A-LR, ensemble: r1i1p1f1, grid: gr, esgf_version: v20180803} + - {institute: KIOST, dataset: KIOST-ESM, ensemble: r1i1p1f1, grid: gr1, esgf_version: v20210601} + - {institute: MIROC, dataset: MIROC-ES2L, ensemble: r1i1p1f2, grid: gn, esgf_version: v20191129} + - {institute: MIROC, dataset: MIROC6, ensemble: r1i1p1f1, grid: gn, esgf_version: v20191016} + - {institute: MPI-M, dataset: MPI-ESM1-2-HR, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190710} + - {institute: MPI-M, dataset: MPI-ESM1-2-LR, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190710} + - {institute: MRI, dataset: MRI-ESM2-0, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190603} + - {institute: NASA-GISS, dataset: GISS-E2-2-G, ensemble: r1i1p1f1, grid: gn, esgf_version: v20191120, start_year: 1970} + - {institute: NCAR, dataset: CESM2-FV2, ensemble: r1i1p1f1, grid: gn, esgf_version: v20191120} + - {institute: NCAR, dataset: CESM2-WACCM-FV2, ensemble: r1i1p1f1, grid: gn, esgf_version: v20191120} + - {institute: NCAR, dataset: CESM2, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190308} + - {institute: NCC, dataset: NorCPM1, ensemble: r1i1p1f1, grid: gn, esgf_version: v20200724} + - {institute: NCC, dataset: NorESM2-LM, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190815} + - {institute: NCC, dataset: NorESM2-MM, ensemble: r1i1p1f1, grid: gn, esgf_version: v20191108} + - {institute: NOAA-GFDL, dataset: GFDL-CM4, ensemble: r1i1p1f1, grid: gr1, esgf_version: v20180701} + - {institute: NOAA-GFDL, dataset: GFDL-CM4, ensemble: r1i1p1f1, grid: gr2, esgf_version: v20180701} + - {institute: NOAA-GFDL, dataset: GFDL-ESM4, ensemble: r1i1p1f1, grid: gr1, esgf_version: v20190726} + - {institute: NUIST, dataset: NESM3, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190812} + - {institute: SNU, dataset: SAM0-UNICON, ensemble: r1i1p1f1, grid: gn, esgf_version: v20190323} + +datasets_ERA5: &datasets_era5 + - {project: OBS6, dataset: ERA5, type: reanaly, frequency: day, latestversion: v1, mip: day, tier: 3, start_year: 1958, end_year: 2014} + +datasets_EOBS: &datasets_EOBS + - {dataset: E-OBS, project: OBS, version: v29.0e-0.25, type: ground, tier: 2, start_year: 1958, end_year: 2014, mip: day} + +preprocessors: + weathertype_preproc: + extract_time: &time + start_year: 1950 + start_month: 1 + start_day: 1 + end_year: 2014 + end_month: 12 + end_day: 31 + regrid: + # This is the region for which the Lamb weathertypes will be calculated. + # The objective classification scheme of Jones et al. 1993 is used, + # these grid points are the same ones as used in the paper. + target_grid: &grid + start_longitude: -5 + end_longitude: 25 + step_longitude: 5 + start_latitude: 35 + end_latitude: 55 + step_latitude: 5 + scheme: linear + + weathertype_preproc_pr: + extract_time: + *time + extract_region: + # This is the region for which precipitation data can be used to + # calculate correlations between weathertypes and this said region + # to further be able to combine them into 'simplified' weathertypes. + # Combined weathertypes correspond to a specific precipitation + # pattern over this region. You can change this to a region of your liking + # to see which types of weathertypes are associated with which precipitation + # patterns. + start_longitude: 9.5 + end_longitude: 17.25 + start_latitude: 46.25 + end_latitude: 49 + + mean_preproc: + extract_time: + *time + extract_region: ® + start_longitude: -15 + end_longitude: 35 + start_latitude: 25 + end_latitude: 65 + +diagnostics: + weathertyping: + description: calculate weathertypes and plot means + variables: + era5_msl_wt: &era5_msl_wt + project: OBS6 + dataset: ERA51 + type: reanaly + frequency: day + mip: day + latestversion: v1 + tier: 3 + start_year: 1958 + end_year: 2014 + short_name: psl + preprocessor: weathertype_preproc + additional_datasets: + *datasets_era5 + + era5_pr_wt: &era5_pr_wt + project: OBS6 + dataset: ERA51 + type: reanaly + frequency: day + mip: day + latestversion: v1 + tier: 3 + start_year: 1958 + end_year: 2014 + short_name: pr + preprocessor: weathertype_preproc_pr + additional_datasets: + *datasets_era5 + + eobs_pr_wt: &eobs_pr_wt + project: OBS + dataset: E-OBS + type: ground + mip: day + version: v29.0e-0.25 + tier: 2 + start_year: 1958 + end_year: 2014 + short_name: pr + preprocessor: weathertype_preproc_pr + additional_datasets: + *datasets_EOBS + + era5_msl_mean: &era5_msl + project: OBS6 + dataset: ERA5 + type: reanaly + frequency: day + mip: day + latestversion: v1 + tier: 3 + start_year: 1958 + end_year: 2014 + short_name: psl + preprocessor: mean_preproc + additional_datasets: + *datasets_era5 + + era5_tp_mean: &era5_prcp + project: OBS6 + dataset: ERA5 + type: reanaly + frequency: day + mip: day + latestversion: v1 + tier: 3 + start_year: 1958 + end_year: 2014 + short_name: pr + preprocessor: mean_preproc + additional_datasets: + *datasets_era5 + + era5_tas_mean: &era5_temp + project: OBS6 + dataset: ERA5 + type: reanaly + frequency: day + mip: day + latestversion: v1 + tier: 3 + start_year: 1958 + end_year: 2014 + short_name: tas + preprocessor: mean_preproc + additional_datasets: + *datasets_era5 + + cmip6_historical_psl_day_wt: &cmip6_historical_wt + project: CMIP6 + activity: CMIP + exp: historical + mip: day + short_name: psl + preprocessor: weathertype_preproc + start_year: 1950 + end_year: 2014 + additional_datasets: + *datasets_psl_cmip6 + + cmip6_historical_psl_day_mean: &cmip6_historical_psl + project: CMIP6 + activity: CMIP + exp: historical + mip: day + short_name: psl + preprocessor: mean_preproc + start_year: 1950 + end_year: 2014 + additional_datasets: + *datasets_psl_cmip6 + + cmip6_historical_prcp_day_mean: &cmip6_historical_prcp + project: CMIP6 + activity: CMIP + exp: historical + mip: day + short_name: pr + preprocessor: mean_preproc + start_year: 1950 + end_year: 2014 + additional_datasets: + *datasets_psl_cmip6 + + cmip6_historical_temp_day_mean: &cmip6_historical_ta + project: CMIP6 + activity: CMIP + exp: historical + mip: day + short_name: tas + preprocessor: mean_preproc + start_year: 1950 + end_year: 2014 + additional_datasets: + *datasets_psl_cmip6 + + scripts: + weathertyping: + correlation_threshold: 0.9 + rmse_threshold: 0.002 + automatic_slwt: true + predefined_slwt: + 1: + - 1 + - 2 + - 19 + 2: + - 3 + - 4 + - 22 + - 21 + 3: + - 5 + - 6 + - 15 + - 16 + 4: + - 7 + - 8 + - 11 + - 18 + 5: + - 9 + - 17 + 6: + - 10 + - 20 + 7: + - 12 + - 13 + - 14 + 8: + - 23 + - 24 + 9: + - 25 + - 26 + 0: + - 27 + plotting: true + script: weathertyping/weathertyping.py diff --git a/esmvaltool/references/jones93ijc.bibtex b/esmvaltool/references/jones93ijc.bibtex new file mode 100644 index 0000000000..c3f0ee590d --- /dev/null +++ b/esmvaltool/references/jones93ijc.bibtex @@ -0,0 +1,14 @@ +@article{jones93ijc, + author = {Jones, P. D. and Hulme, M. and Briffa, K. R.}, + title = {A comparison of Lamb circulation types with an objective classification scheme}, + journal = {International Journal of Climatology}, + volume = {13}, + number = {6}, + pages = {655-663}, + keywords = {Weather types, Temperature, Precipitation, British, Isles}, + doi = {https://doi.org/10.1002/joc.3370130606}, + url = {https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/joc.3370130606}, + eprint = {https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/joc.3370130606}, + abstract = {Abstract An objective scheme, initially developed by Jenkinson and Collison, is used to classify daily circulation types over the British Isles, along the lines of the subjective method devised by Lamb. The scheme uses daily grid-point mean sea-level pressure data for the region. The results of the analysis over the period 1881–1989 are compared with ‘true’ Lamb weather types. The frequencies of objectively developed types are highly correlated with traditional Lamb types, especially so for synoptic (cyclonic and anticyclonic) types, although still good for wind directional types. Comparison of the two classification schemes reveals negligible differences between the correlations of the counts of circulation types and regional temperature and rainfall. The major difference between the two classification schemes is that the decline of the westerlies since 1940 is less evident with the objective scheme.}, + year = {1993} +} diff --git a/esmvaltool/references/maraun21jgr.bibtex b/esmvaltool/references/maraun21jgr.bibtex new file mode 100644 index 0000000000..d160c9e78f --- /dev/null +++ b/esmvaltool/references/maraun21jgr.bibtex @@ -0,0 +1,15 @@ +@article{maraun21jgr, + author = {Maraun, Douglas and Truhetz, Heimo and Schaffer, Armin}, + title = {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 = {Journal of Geophysical Research: Atmospheres}, + volume = {126}, + number = {6}, + pages = {e2020JD032824}, + keywords = {Austria, bias adjustment, climate model evaluation, large-scale circulation errors, regional climate projections}, + doi = {https://doi.org/10.1029/2020JD032824}, + url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020JD032824}, + eprint = {https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020JD032824}, + note = {e2020JD032824 2020JD032824}, + abstract = {Abstract The Austrian regional climate projections are based on an ensemble of bias adjusted regional climate model simulations. Bias adjustment (BA) improves the usability of climate projections for impact studies, but cannot mitigate fundamental model errors. This argument holds in particular for biases in the temporal dependence, which is strongly influenced by the large-scale circulation. Global climate models (GCMs), underlying regional climate projections, suffer from substantial circulation errors. We therefore, conduct a process-based evaluation of the Austrian climate projections focusing on large-scale circulation errors, their regional imprints and the potential for BA. First, we define nine synoptic weather types and assess how well the considered climate models represent their occurrence and persistence. Second, we assess biases in the overall dry and hot day probabilities, as well as conditional on synoptic weather type occurrence; and biases in the duration of dry and hot spells. Third, we investigate how these biases depend on biases in the occurrence and persistence of relevant synoptic weather types. And fourth, we study how much an overall BA improves these biases. Many GCMs misrepresent the occurrence and persistence of relevant synoptic weather types. These biases have a clear imprint on biases in dry and hot day occurrence and spell durations. BA in many cases helps to greatly reduce biases even in the presence of circulation biases, but may in some cases amplify conditional biases. Persistence biases are especially relevant for the representation of meteorological drought. Biases in the duration of dry spells cannot fully be mitigated by BA.}, + year = {2021} +}