diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 9eec648279..4e2698454f 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -6,7 +6,6 @@ on: push: branches: - main - - fix_recipe_filler_bkwds_incompatibility schedule: - cron: '0 0 * * *' diff --git a/CITATION.cff b/CITATION.cff index 5c253e3bb5..cd621538b7 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -175,6 +175,11 @@ authors: family-names: Hagemann given-names: Stefan orcid: "https://orcid.org/0000-0001-5444-2945" + - + affiliation: "University of Canterbury, New Zealand" + family-names: Hardacre + given-names: Catherine + orcid: "https://orcid.org/0000-0001-9093-4656" - affiliation: "ISAC-CNR, Italy" name-particle: von diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_DJF.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_DJF.png new file mode 100644 index 0000000000..ccf5c4b1c8 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_DJF.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_JJA.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_JJA.png new file mode 100644 index 0000000000..38b09c52f9 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_JJA.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_MAM.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_MAM.png new file mode 100644 index 0000000000..5e7ca417a1 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_MAM.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_SON.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_SON.png new file mode 100644 index 0000000000..ffb666eea2 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_SON.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_scatter.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_scatter.png new file mode 100644 index 0000000000..67c64a3fa7 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1988_2008_scatter.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_DJF.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_DJF.png new file mode 100644 index 0000000000..f653410e07 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_DJF.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_JJA.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_JJA.png new file mode 100644 index 0000000000..6474acc856 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_JJA.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_MAM.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_MAM.png new file mode 100644 index 0000000000..2e2fda4cca Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_MAM.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_SON.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_SON.png new file mode 100644 index 0000000000..ce4d3fac08 Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_SON.png differ diff --git a/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_scatter.png b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_scatter.png new file mode 100644 index 0000000000..9226bca81a Binary files /dev/null and b/doc/sphinx/source/recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_scatter.png differ diff --git a/doc/sphinx/source/recipes/index.rst b/doc/sphinx/source/recipes/index.rst index edcc48977a..0f0ce7667d 100644 --- a/doc/sphinx/source/recipes/index.rst +++ b/doc/sphinx/source/recipes/index.rst @@ -62,6 +62,7 @@ Atmosphere recipe_thermodyn_diagtool recipe_validation recipe_radiation_budget + recipe_aod_aeronet_assess Climate metrics ^^^^^^^^^^^^^^^ diff --git a/doc/sphinx/source/recipes/recipe_aod_aeronet_assess.rst b/doc/sphinx/source/recipes/recipe_aod_aeronet_assess.rst new file mode 100644 index 0000000000..fec1bed761 --- /dev/null +++ b/doc/sphinx/source/recipes/recipe_aod_aeronet_assess.rst @@ -0,0 +1,161 @@ +.. _recipe_aod_aeronet_assess: + +AOD AeroNET Assess +================== + +Overview +-------- + +This diagnostic evaluates model aerosol optical depth (AOD) against ground +based observations from the AeroNET measurement network. Monthly mean AOD +data is downloaded from the AeroNET website and formatted (CMORized) using the +AERONET downloader and formatter within ESMValTool. + +Multiannual seasonal means are calculated from the model output and compared +with a multiannual seasonal mean climatology generated from AeroNET +observational data. At each AeroNET station the data are screened for validity +according to the following default criteria: + + * 1. Monthly means must be generated from at least one AOD observation in that + month. + + * 2. Seasonal means for DJF, MAM, JJA and SON must be calculated from three + monthly means, i.e. a monthly mean from December January and Feburary. + + * 3. For a given year to be valid, there must be a seasonal mean for each climate + season i.e. DJF, MAM, JJA and SON. + + * 4. For a multiannual seasonal means there must be at least five seasonaal means + over the time range of interest. + +NOTE: The code is designed to be flexible and the default criteria can be +changed according to the user's requirements (see the user settings below). + +The evaluation is visualised by plotting model output as 2D filled contours and +overlaying AeroNET observations at model grid cells co-located with the AeroNET +measurement stations. Statistical data (root mean square error) is generated +using AeroNET observations at model grid cells co-located with the AeroNET +measurement stations. + +Available recipes and diagnostics +--------------------------------- + +Recipes are stored in esmvaltool/recipes/ + + * recipe_aod_aeronet_assess.yml + +Diagnostics are stored in esmvaltool/diag_scripts/aerosols/ + + * aod_aeronet_assess.py: Plot the AOD evaluation. + * aero_utils.py: Utility functions commonly used by aerosol assessment routines. + + +User settings in recipe +----------------------- + +#. Script aod_aeronet_assess.py + + *Required settings for script* + + * wavel: The wavelength of interest for the evaluation, currently set up for 440nm + * min_days_per_mon: The minimum number of days used to calculate the AOD monthly mean + * min_mon_per_seas: The minimum number of seasons used to calculate each + seasonal mean. This must be between 1 and 3. + * min_seas_per_year: The minimum number of seasonal means in each year. This + must be between 1 and 4. + * min_seas_per_clim: The minimum number of seasonal means used to calculate + the multiannual seasonal mean. This must be btween 1 and the number of years + of available AeroNET data. + + *Optional settings for script* + + * None + + *Required settings for variables* + + * None + + *Optional settings for variables* + + * None + + *Required settings for preprocessor* + + * None + + *Optional settings for preprocessor* + + * None + + *Color tables* + + * brewer_Spectral_11 + + +Variables +--------- + +* od440aer (atmos, monthly mean, longitude latitude time) + + +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. + +* The AeroNET data is downloaded from the AeroNET website using the downloader: + + .. code-block:: yaml + + $ esmvaltool data download AERONET. + +* The AeroNET data is formatteed (CMORized) using the formatter: + + .. code-block:: yaml + + $ esmvaltool data format AERONET. + + + +References +---------- +* Holben B.N., T.F.Eck, I.Slutsker, D.Tanre, J.P.Buis, A.Setzer, E.Vermote, J.A.Reagan, Y.Kaufman, T.Nakajima, F.Lavenu, I.Jankowiak, and A.Smirnov, 1998: AERONET - A federated instrument network and data archive for aerosol characterization, Rem. Sens. Environ., 66, 1-16. + +* Holben, B.N., D.Tanre, A.Smirnov, T.F.Eck, I.Slutsker, N.Abuhassan, W.W.Newcomb, J.Schafer, B.Chatenet, F.Lavenue, Y.J.Kaufman, J.Vande Castle, A.Setzer, B.Markham, D.Clark, R.Frouin, R.Halthore, A.Karnieli, N.T.O'Neill, C.Pietras, R.T.Pinker, K.Voss, and G.Zibordi, 2001: An emerging ground-based aerosol climatology: Aerosol Optical Depth from AERONET, J. Geophys. Res., 106, 12 067-12 097. + +* Mulcahy, J. P., Johnson, C., Jones, C. G., Povey, A. C., Scott, C. E., Sellar, A., Turnock, S. T., Woodhouse, M. T., Abraham, N. L., Andrews, M. B., Bellouin, N., Browse, J., Carslaw, K. S., Dalvi, M., Folberth, G. A., Glover, M., Grosvenor, D. P., Hardacre, C., Hill, R., Johnson, B., Jones, A., Kipling, Z., Mann, G., Mollard, J., O’Connor, F. M., Palmiéri, J., Reddington, C., Rumbold, S. T., Richardson, M., Schitgens, N. A. J., Stier, P., Stringer, M., Tang, Y., Walton, J., Woodward, S., and Yool. A.: Description and evaluation of aerosol in UKESM1 and HadGEM3-GC3.1 CMIP6 historical simulations, Geosci. Model Dev., 13, 6383–6423, 2020 + +Example plots +------------- + +.. _fig_aod_aeronet_assess_1: +.. figure:: /recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_DJF.png + :align: center + + Evaluation of AOD at 440 nm from UKESM1 historical ensemble member r1i1p1f2 against the AeroNET climatology from ground-based observations for Dec-Jan-Feb. The multiannual seasonal mean is calculated for the model data for the period 1994-2014. The model output is overlaid with the observational climatology. + +.. _fig_aod_aeronet_assess_2: +.. figure:: /recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_MAM.png + :align: center + + Evaluation of AOD at 440 nm from UKESM1 historical ensemble member r1i1p1f2 against the AeroNET climatology from ground-based observations for Mar_Apr_May. The multiannual seasonal mean is calculated for the model data for the period 1994-2014. The model output is overlaid with the observational climatology. + +.. _fig_aod_aeronet_assess_3: +.. figure:: /recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_JJA.png + :align: center + + Evaluation of AOD at 440 nm from UKESM1 historical ensemble member r1i1p1f2 against the AeroNET climatology from ground-based observations for Jun-Jul-Aug. The multiannual seasonal mean is calculated for the model data for the period 1994-2014. The model output is overlaid with the observational climatology. + +.. _fig_aod_aeronet_assess_4: +.. figure:: /recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_SON.png + :align: center + + Evaluation of AOD at 440 nm from UKESM1 historical ensemble member r1i1p1f2 against the AeroNET climatology from ground-based observations for Sep-Oct-Nov. The multiannual seasonal mean is calculated for the model data for the period 1994-2014. The model output is overlaid with the observational climatology. + +.. _fig_aod_aeronet_assess_5: +.. figure:: /recipes/figures/aod_aeronet_assess/UKESM1-0-LL_CMIP_AERmon_historical_od440aer_gn_1994_2014_scatter.png + :align: center + + Evaluation of AOD at 440 nm from UKESM1 historical ensemble member r1i1p1f2 against the AeroNET climatology from ground-based observations for Dec-Jan-Feb, Mar_Apr_May, Jun-Jul-Aug and Sep-Oct-Nov. The multiannual seasonal mean is calculated for the model data for the period 1994-2014. diff --git a/esmvaltool/config-references.yml b/esmvaltool/config-references.yml index 285ace740e..b5f43bc911 100644 --- a/esmvaltool/config-references.yml +++ b/esmvaltool/config-references.yml @@ -252,6 +252,10 @@ authors: name: Hansson, Ulf institute: SMHI, Sweden orcid: + hardacre_catherine: + name: Hardacre, Catherine + institute: University of Canterbury, New Zealand + orcid: https://orcid.org/0000-0001-9093-4656 hassler_birgit: name: Hassler, Birgit institute: DLR, Germany diff --git a/esmvaltool/diag_scripts/aerosols/aero_utils.py b/esmvaltool/diag_scripts/aerosols/aero_utils.py new file mode 100644 index 0000000000..623a85f2aa --- /dev/null +++ b/esmvaltool/diag_scripts/aerosols/aero_utils.py @@ -0,0 +1,193 @@ +"""Part of the ESMValTool Aerosol diagnostics. + +This module contains utility functions commonly used by aerosol +assessment routines. +""" + +import iris +import numpy as np + + +class AeroAnsError(Exception): + + """Exception class for errors raised when model data is checked in the + extract_pt module. + """ + + +def add_bounds(cube): + """Add bounds to a cube's latitude and longitude coordinates. + + Parameters + ---------- + cube : Iris cube + Iris cube with latitude and longitude coordinates. + + Returns + ------- + cube : Iris cube. + Iris cube with bounds added to the latitude and longitude coordinates. + """ + + if not cube.coord('latitude').has_bounds(): + cube.coord('latitude').guess_bounds() + if not cube.coord('longitude').has_bounds(): + cube.coord('longitude').guess_bounds() + + return cube + + +def extract_pt(icube, pt_lat, pt_lon, height=None, level=None, nearest=False): + """Extracts given location(s) (3-D) from a cube. + + Method + ------ + Uses Iris module Analysis.Interpolate to extract values, + initially based on horizontal coordinates, and then based on + height, if specified. + + If height ('altitude') is requested, checks if cube heights + include orography, i.e. HybridHeights have been derived. + + Parameters + ---------- + icube : Iris cube + pt_lat, pt_lon : Float or list/array of floats. Latitude and longitude + coordinates of desired points. + args: + height : Float or list/array of floats. Altitude (above geoid) of + point. Initialized to None. + level : Integer . Model level or pseudo level or tile number. + Initialized to None, meaning that all available levels in + the cube are used. + nearest : Boolean. Specify whether to use 'nearest neighbour', instead + of 'linear' method while extracting data. Default is False. + + Returns + ------- + data_out : List + List of single point values, corresponding to each point specified. + + Raises + ------ + AeroAnsError : If the number of latitude and longitude points are + mismatched. OR if both and level and height are passed as args. + OR if the cube contains a time coordinate. OR if a pseudo level + coordinate is requested, but not present in the cube. OR if the numbers + of latitude/longitude and height points are mismatched. OR if height + is requested but the cube does not contain an altitude coordinate. + """ + + # Check that input data is a (single) cube + if not isinstance(icube, iris.cube.Cube): + raise AeroAnsError('Extract_pt:First argument must be a single cube') + + # Check if the cube contains a time dimension, which is + # currently unsupported. + if icube.coords()[0].name() == 'time': + raise AeroAnsError( + 'Extract_pt:Cannot handle time dimension at present') + + # Check that equal number of lat/lon pairs are passed in point coordinates. + # Convert arguments to lists for easier processing if necessary. + pt_lat1 = [] + pt_lon1 = [] + + if not isinstance(pt_lat, list): + pt_lat1.append(pt_lat) + pt_lon1.append(pt_lon) + + else: + for n_lat in np.arange(len(pt_lat)): + pt_lat1.append(pt_lat[n_lat]) + pt_lon1.append(pt_lon[n_lat]) + + if len(pt_lat1) != len(pt_lon1): + raise AeroAnsError('Extract_pt:Mismatch in number of lat/long values') + + # Check that both level and height haven't been requested. + if level is not None and height is not None: + raise AeroAnsError('Extract_pt: Both Level and Height requested') + + # Check that the cube has a level coordinate if level has been requested. + if level is not None and not icube.coord( + 'model_level_number') and not icube.coord('pseudo_level'): + raise AeroAnsError('Extract_pt:Level requested, but not found in cube') + + # Check that the number of height points is equal to the number of + # lat/lon pairs. Convert the argument to a list for easier + # processing if necessary. + if height is not None: + pt_hgt = [] + +# if isinstance(height, list): +# pt_hgt.extend(height) +# else: +# pt_hgt.append(height) + pt_hgt.extend(height) if \ + isinstance(height, list) else \ + pt_hgt.append(height) + + if len(pt_lat1) != len(pt_hgt): + raise AeroAnsError( + 'Extract_pt:Mismatch in number of points for lat/long/height') + + # Check that heights have been merged with orography. + if not icube.coords('altitude'): + raise AeroAnsError( + 'Extract_pt:Height requested but input data does not contain \ + "Altitude" coordinate') + + # Store the min and max altitudes from cube data so that user + # cannot request points located below/ above that. + # Will extract =min/max if beyond limits. + hgt_min = icube.coord('altitude').points.min() + hgt_max = icube.coord('altitude').points.max() + + # ---------- Finished checks -- begin processing ------------------------- + + # If level specified, extract slice first + if level is not None: + + try: + icube = icube.extract( + iris.Constraint(model_level_number=level)) + + except Exception: + print('Model level number not available. Use pseudo level') + + else: + icube = icube.extract( + iris.Constraint(pseudo_level=level)) + + # Extract values for specified points lat/lon + # NOTE: Does not seem to handle multiple points if 3-D + data_out = [] + + # Set lat/lon coordinates for model grid cell interpolation + for n_lat1 in np.arange(len(pt_lat1)): + latlon_coords = [('latitude', pt_lat1[n_lat1]), + ('longitude', pt_lon1[n_lat1])] + + if nearest: + tcube = icube.interpolate(latlon_coords, iris.analysis.Nearest()) + else: + tcube = icube.interpolate(latlon_coords, iris.analysis.Linear()) + + # If height specified, interpolate to requested height + if height is not None: + + # Set vertical coordinates for model grid cell interpolation + point = max(pt_hgt[n_lat1], hgt_min) + point = min(pt_hgt[n_lat1], hgt_max) + hgt_coords = [('altitude', point)] + + if nearest: + tcube = tcube.interpolate(hgt_coords, iris.analysis.Nearest()) + else: + tcube = tcube.interpolate(hgt_coords, iris.analysis.Linear()) + + # Append processed data point + data_out.append(tcube.data) + + return data_out diff --git a/esmvaltool/diag_scripts/aerosols/aod_aeronet_assess.py b/esmvaltool/diag_scripts/aerosols/aod_aeronet_assess.py new file mode 100644 index 0000000000..27ab6b2714 --- /dev/null +++ b/esmvaltool/diag_scripts/aerosols/aod_aeronet_assess.py @@ -0,0 +1,448 @@ +"""Implement the AOD climatology metric from ground-based + AeroNet observations. +""" +import logging +import os + +import iris +import iris.plot as iplt +import matplotlib.cm as mpl_cm +import matplotlib.lines as mlines +import matplotlib.pyplot as plt +import numpy as np +import scipy +from aero_utils import add_bounds, extract_pt +from matplotlib import colors, gridspec +from numpy import ma + +from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic +from esmvaltool.diag_scripts.shared._base import get_plot_filename + +logger = logging.getLogger(os.path.basename(__file__)) +fontsizedict = {"title": 25, "axis": 20, "legend": 18, "ticklabel": 18} + + +def get_provenance_record(filenames): + """Return a provenance record describing the metric. + + Parameters + ---------- + filenames : List of strings + The filenames containing the data used to create the metric. + + Returns + ------- + dictionary + The provenance record describing the metric. + """ + record = { + "ancestors": filenames, + } + + return record + + +def plot_aod_mod_obs(md_data, obs_data, aeronet_obs_cube, plot_dict): + """Plot AOD contour overlaid with Aeronet climatology. + + Parameters + ---------- + md_data : Iris cube + Model AOD as a cube with latitude and longitude coordinates. + obs_data : List. + Observations of AOD from each AeroNET station. + aeronet_obs_cube : Iris cube. + Holds information about Aeronet measurement stations including + station names, station latitude and station longitude. + plot_dict : Dictionary. + Holds plotting settings. + """ + # Plot model data + cf_plot = iplt.contourf(md_data, + plot_dict["Levels"], + colors=plot_dict["Colours"], + extend="max") + + # Latitude and longitude of stations. + anet_aod_lats = aeronet_obs_cube.coord("latitude").points + anet_aod_lons = ((aeronet_obs_cube.coord("longitude").points + 180) % 360 - + 180) + + # Loop over stations + for istn, stn_data in enumerate(obs_data): + if ma.is_masked(stn_data): + continue + + # Find position of the observed AOD on the colorscale. + # np.searchsorted returns index at which inserting new value will + # maintain a sorted array. We use the color to the left of index. + cid = np.searchsorted(plot_dict["Levels"], stn_data) + cid = max(0, cid - 1) # filter out zero and max when seeking 'left' + cid = min(len(plot_dict["Colours"]) - 1, cid) + pcol = plot_dict["Colours"][cid] + + # Overlay contourf with observations + plt.plot( + anet_aod_lons[istn], + anet_aod_lats[istn], + color=pcol, + marker="o", + markeredgecolor="k", + markeredgewidth=2, + markersize=9, + ) + + # Decorate the plot + plt.title(plot_dict["Title"], size=24) + colbar = plt.colorbar(cf_plot, orientation="horizontal") + colbar.set_ticks(plot_dict["Levels"]) + colbar.set_ticklabels(plot_dict["tick_labels"]) + plt.gca().coastlines(color="#525252") + + # Statistics on plot + plt.figtext( + 0.12, + 0.27, + (f'''Global mean AOD={plot_dict["Mean_aod"]:.3f}; RMSE=''' + f'''{plot_dict["RMS_aod"]:.3f}; Stn mean: md=''' + f'''{plot_dict["Stn_mn_md"]:.3f}; obs=''' + f'''{plot_dict["Stn_mn_obs"]:.3f}'''), + size=16, + ) + + +def aod_analyse(model_data, aeronet_obs_cube, clim_seas, wavel): + """Evaluates AOD vs Aeronet, generates plots and returns evaluation + metrics. + + Parameters + ---------- + model_data : Iris Cube. + Contains model output of AOD with coordinates; time, latitude and + longitude. + aeronet_obs_cube : Iris Cube. + Contains information about Aeronet measurement stations including + station names, station latitude and station longitude. + clim_seas : List. + Strings to denote climate seasons ["DJF", "MAM", "JJA", "SON"] + wavel : String. + AOD wavelength, default = 440nm - translates to pseudo-level. + + Returns + ------- + figures : List. + Contains figure instances for the seasonal contour plots overlaid with + observations of AOD from AeroNET. + fig_scatter : Figure object. + The scatter plot comparing modelled and observed AOD at 440nm. + """ + # Convert wave length nm -> um + wv_mi = str(float(wavel) / 1000.0) + + # Get model run id + if "parent_source_id" in model_data.attributes: + model_id = model_data.attributes["parent_source_id"] + else: + model_id = "Multi-Model-Mean" + + # Add bounds for lat and lon if not present + model_data = add_bounds(model_data) + + # Co-locate model grid points with measurement sites --func from aero_utils + anet_aod_lats = aeronet_obs_cube.coord("latitude").points.tolist() + anet_aod_lons = aeronet_obs_cube.coord("longitude").points.tolist() + aod_at_anet = extract_pt(model_data, anet_aod_lats, anet_aod_lons) + + # Set up seasonal contour plots + figures = [] + + clevs = [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 2.0] + clabs = [ + "0.0", "", "0.1", "", "0.2", "", "0.3", "", "0.4", "", "0.5", "2.0" + ] + cmapr = mpl_cm.get_cmap("brewer_Spectral_11") + cmap = colors.ListedColormap(cmapr(range(cmapr.N))[-1::-1]) + colours = cmap.colors + + # Set up the figure for scatter plotting + fig_scatter = plt.figure(figsize=(10, 10)) + gs_scatter = gridspec.GridSpec(ncols=1, nrows=1) + ax_scatter = fig_scatter.add_subplot(gs_scatter[0, 0]) + col_scatter = ["#081d58", "#41ab5d", "#fe9929", "#7f0000"] + leg_scatter = [] + + # Loop over seasons + for season in aeronet_obs_cube.slices_over("clim_season"): + + # Match Aeronet obs season with model season number + model_sn = [c.lower() for c in clim_seas + ].index(season.coord("clim_season").points[0]) + model_season = model_data[model_sn] + + logger.info('Analysing AOD for %s: %s', {model_id}, + {clim_seas[model_sn]}) + + # Generate statistics required - area-weighted mean + grid_areas = iris.analysis.cartography.area_weights(model_season) + global_mean = model_season.collapsed( + ["latitude", "longitude"], + iris.analysis.MEAN, + weights=grid_areas, + ) + + # Extract model and obs data for season number (model_sn) + seas_anet_obs = season.data + seas_anet_md = np.array([x[model_sn] for x in aod_at_anet]) + + # Match model data with valid obs data + valid_indices = ma.where(seas_anet_obs) + valid_obs = seas_anet_obs[valid_indices] + valid_md = seas_anet_md[valid_indices] + + # Model - obs statistics (diff, model mean and RMS, r2) + diff = valid_md - valid_obs + stn_mn_obs = np.mean(valid_obs) + stn_mn_md = np.mean(valid_md) + rms_aod = np.sqrt(np.mean(diff**2)) + linreg = scipy.stats.linregress(valid_obs, valid_md) + + # Plot scatter of co-located model and obs data + ax_scatter.scatter(valid_obs, valid_md, color=col_scatter[model_sn]) + + # Legend + label = f"{clim_seas[model_sn]} = {linreg.rvalue**2:.2f}" + leg_scatter.append( + mlines.Line2D( + [0], + [0], + marker="o", + color="w", + label=label, + markersize=15, + markerfacecolor=col_scatter[model_sn], + )) + + # Plot contours overlaid with obs for this run and season + fig_cf = plt.figure(figsize=(11, 8), dpi=300) + + n_stn = str(len(valid_obs)) + title = ("\nTotal Aerosol Optical Depth at " + wv_mi + " microns" + + "\n" + model_id + ", " + clim_seas[model_sn] + + ", N stations=" + n_stn) + + # Plot dictionary + plot_dict = { + "Mean_aod": global_mean.data, + "Stn_mn_obs": stn_mn_obs, + "Stn_mn_md": stn_mn_md, + "RMS_aod": rms_aod, + "Levels": clevs, + "Colours": colours, + "tick_labels": clabs, + "Title": title, + "Season": clim_seas[model_sn], + } + plot_aod_mod_obs(model_season, seas_anet_obs, aeronet_obs_cube, + plot_dict) + + figures.append(fig_cf) + + # Decorate the scatter plot + line = mlines.Line2D([0, 1], [0, 1], color="#696969") + transform = ax_scatter.transAxes + line.set_transform(transform) + ax_scatter.add_line(line) + + ax_scatter.set( + xlim=(0, 1), + xticks=np.linspace(0.0, 1.0, num=6), + ylim=(0, 1), + yticks=np.linspace(0.0, 1.0, num=6), + ) + ax_scatter.set_xlabel("AeroNET AOD", fontsize=fontsizedict["axis"]) + ax_scatter.set_ylabel(model_id + " AOD", fontsize=fontsizedict["axis"]) + + ax_scatter.tick_params(axis="both", + which="major", + labelsize=fontsizedict["ticklabel"]) + + ax_scatter.set_title( + "Model vs obs: Total Aerosol Optical Depth \n at " + wv_mi + + " microns", + fontsize=fontsizedict["title"], + ) + + ax_scatter.legend( + handles=leg_scatter, + loc="lower right", + title="Seasonal R2", + title_fontsize=fontsizedict["legend"], + fontsize=fontsizedict["legend"], + ) + + return figures, fig_scatter + + +def preprocess_aod_obs_dataset(obs_dataset): + """Calculate a multiannual seasonal mean AOD climatology. + + Observational AOD timeseries data from AeroNET are used to generate a + multiannual seasonal mean climatology for each AeroNET station. The + user sets thresholds (or uses the default settings) to specify the + amount of valid data required for the climatology. At this stage + ESMValTool preprocessors are unsuitable for pre-processing the AeroNET + AOD observations because of the bespoke nature and application of the + filtering thresholds. + + Parameters + ---------- + obs_dataset : ESMValTool dictionary. Holds meta data for the observational + AOD dataset. + + Returns + ------- + multiannual_seaonal_mean : Iris cube. Preprocessed observational + AOD climatology. + """ + obs_cube = iris.load_cube(obs_dataset[0]["filename"]) + + # Set up thresholds for generating the multi annual seasonal mean + min_days_per_mon = 1 + min_mon_per_seas = 3 + min_seas_per_year = 4 + min_seas_per_clim = 5 + + # Add the clim_season and season_year coordinates. + iris.coord_categorisation.add_year(obs_cube, 'time', name='year') + + iris.coord_categorisation.add_season(obs_cube, 'time', name='clim_season') + + iris.coord_categorisation.add_season_year(obs_cube, + 'time', + name='season_year') + + # Copy obs cube and mask all months with fewer + # "Number of days" than given threshold. + num_days_var = obs_cube.ancillary_variable("Number of days") + masked_months_obs_cube = obs_cube.copy(data=ma.masked_where( + num_days_var.data < min_days_per_mon, obs_cube.data)) + + # Aggregate (mean) by season. + # The number of unmasked months per season is counted, + # and where there are fewer unmasked months than the + # given threshold, the computed mean is masked. + annual_seasonal_mean = masked_months_obs_cube.aggregated_by( + ['clim_season', 'season_year'], + iris.analysis.MEAN, + ) + annual_seasonal_count = masked_months_obs_cube.aggregated_by( + ['clim_season', 'season_year'], + iris.analysis.COUNT, + function=lambda values: ~ma.getmask(values), + ) + annual_seasonal_mean.data = ma.masked_where( + annual_seasonal_count.data < min_mon_per_seas, + annual_seasonal_mean.data, + ) + + # Aggregate (mean) by multi-annual season. + # The number of unmasked seasons per multi-annual season + # is counted, and where there are fewer unmasked seasons + # than the given threshold, the computed multi-annual + # season is masked. + multi_annual_seasonal_mean = annual_seasonal_mean.aggregated_by( + 'clim_season', + iris.analysis.MEAN, + ) + clim_season_agg_count = annual_seasonal_mean.aggregated_by( + 'clim_season', + iris.analysis.COUNT, + function=lambda values: ~ma.getmask(values), + ) + multi_annual_seasonal_mean.data = ma.masked_where( + clim_season_agg_count.data < min_seas_per_clim, + multi_annual_seasonal_mean.data, + ) + year_agg_count = multi_annual_seasonal_mean.aggregated_by( + 'year', + iris.analysis.COUNT, + function=lambda values: ~ma.getmask(values), + ) + + counter = range(len( + multi_annual_seasonal_mean.coord('clim_season').points)) + for iseas in counter: + multi_annual_seasonal_mean.data[iseas, :] = ma.masked_where( + year_agg_count.data[0, :] < min_seas_per_year, + multi_annual_seasonal_mean.data[iseas, :], + ) + + return multi_annual_seasonal_mean + + +def main(config): + """Produce the AOD climatology metric from ground-based AeroNet + observations. + + Parameters + ---------- + wavel : String. + User defined. Default is "440". + config : dict + The ESMValTool configuration. + """ + input_data = config["input_data"] + datasets = group_metadata(input_data.values(), "dataset") + + # Default wavelength + wavel = "440" + + # Produce climatology for observational dataset + obs_dataset = datasets.pop(config["observational_dataset"]) + obs_cube = preprocess_aod_obs_dataset(obs_dataset) + + for model_dataset, group in datasets.items(): + # 'model_dataset' is the name of the model dataset. + # 'group' is a list of dictionaries containing metadata. + logger.info("Processing data for %s", model_dataset) + logger.info(group) + + for attributes in group: + logger.info(attributes["filename"]) + + input_file = attributes["filename"] + provenance_record = get_provenance_record(input_file) + logger.info(provenance_record) + cube = iris.load_cube(input_file) + + # Set up for analysis and plotting + seasons = ["DJF", "MAM", "JJA", "SON"] + + plot_file_prefix = (model_dataset + "_" + attributes["activity"] + + "_" + attributes["mip"] + "_" + + attributes["exp"] + "_" + + attributes["short_name"] + "_" + + str(attributes["start_year"]) + "_" + + str(attributes["end_year"]) + "_") + + # Analysis and plotting for model-obs comparison + figures, fig_scatter = aod_analyse(cube, + obs_cube, + seasons, + wavel=wavel) + + # Save the scatter plot + output_file = plot_file_prefix + "scatter" + output_path = get_plot_filename(output_file, config) + fig_scatter.savefig(output_path) + + # Save the contour plots + for ifig, seas_fig in enumerate(figures): + output_file = plot_file_prefix + seasons[ifig] + output_path = get_plot_filename(output_file, config) + seas_fig.savefig(output_path) + + +if __name__ == "__main__": + with run_diagnostic() as CONFIG: + main(CONFIG) diff --git a/esmvaltool/recipes/recipe_aod_aeronet_assess.yml b/esmvaltool/recipes/recipe_aod_aeronet_assess.yml new file mode 100644 index 0000000000..0fc82ec864 --- /dev/null +++ b/esmvaltool/recipes/recipe_aod_aeronet_assess.yml @@ -0,0 +1,65 @@ +# ESMValTool +# recipe_aod_aeronet_assess.yml +--- +documentation: + description: | + Recipe to plot seasonal maps of global aerosol optical depth (AOD) at 440nm. + + title: Recipe that runs an AOD diagnostic + + authors: + - hogan_emma + - lillis_jon + - hardacre_catherine + + maintainer: + - hogan_emma + - lillis_jon + - hardacre_catherine + + projects: + - esmval + +preprocessors: + ma_season_mean: + regrid: + target_grid: 2.5x2.5 + scheme: nearest + climate_statistics: + operator: mean + period: season + seasons: ['DJF', 'MAM', 'JJA', 'SON'] + multi_model_statistics: + span: overlap + statistics: [mean] + +diagnostics: + od440aer_climatologies: + description: Visualise spatial multi-annual seasonal means AOD at 440nm. + variables: + + od440aer: &var_od440aer + mip: AERmon + short_name: od440aer + start_year: 1994 + end_year: 2014 + additional_datasets: + - {dataset: AERONET, project: OBS6, mip: AERmon, tier: 3, type: atmos, version: 20231021} + + od440aer_season: + <<: *var_od440aer + preprocessor: ma_season_mean + additional_datasets: + - {dataset: UKESM1-0-LL, project: CMIP6, mip: AERmon, exp: historical, ensemble: r1i1p1f2, grid: gn} + - {dataset: HadGEM3-GC31-LL, project: CMIP6, mip: AERmon, exp: historical, ensemble: r1i1p1f3, grid: gn} + - {dataset: EC-Earth3-AerChem, project: CMIP6, mip: AERmon, exp: historical, ensemble: r1i1p1f1, grid: gn} +# - {dataset: NorESM2-LM, project: CMIP6, mip: AERmon, exp: historical, ensemble: r1i1p1f1, grid: gn} + - {dataset: GFDL-ESM4, project: CMIP6, mip: AERmon, exp: historical, ensemble: r1i1p1f1, grid: gr1} + - {dataset: MPI-ESM-1-2-HAM, project: CMIP6, mip: AERmon, exp: historical, ensemble: r1i1p1f1, grid: gn} + + scripts: + aeronet: + script: aerosols/aod_aeronet_assess.py + observational_dataset: AERONET + quickplot: + plot_type: plot