<jupyter_start><jupyter_text>Australian forecast skillThis notebook contains plots showing the skill of quantities averaged over the Natural Resource Management (NRM) super-cluster regions (see https://www.climatechangeinaustralia.gov.au/en/overview/methodology/nrm-regions/). As for the generic skill results, plots are shown for CAFE-f6 and for the CanESM5 and EC-Earth3 CMIP6 DCPP submissions (where possible). 30-year climatological and verification periods are used for both the CAFE-f6 and CMIP6 data. However, because the historical CMIP6 data end in 2014, these periods differ slightly: 1991-2020 for CAFE-f6; 1985-2014 for CanESM5 and EC-Earth3. Note also that the CanESM5 hindcasts are initialised at the end of December every year, while the CAFE-f6 and EC-Earth3 hindcasts are initialised at the beginning of November every year.Significant points (as identified by a non-parametric block bootstrapping procedure) are shown by filled circles for the line plots.Much of this work is motivated by deliverables for the Australian Climate Services (ACS) project. If you wish to see anything added to these results, please contact Dougie.<jupyter_code>import warnings

warnings.filterwarnings("ignore")

import geopandas

import xarray as xr

import matplotlib.pyplot as plt

from src import utils, plot

from notebook_helper import plot_metrics, plot_metric_maps, plot_hindcasts
%load_ext autoreload
%autoreload 2
%load_ext lab_black<jupyter_output><empty_output><jupyter_text>Focus regionsIn the following we'll show some skill results for the following NRM super clusters and for the average over all of Australia. Note, the grid below shows the native CAFE grid.<jupyter_code>fig, ax = plt.subplots(1, 1, figsize=(8, 8))
df = geopandas.read_file("../data/raw/NRM_super_clusters/NRM_super_clusters.shp")
_ = df.plot("label", ax=ax, legend=True, legend_kwds={"loc": "lower left"})

data = xr.open_zarr(
    "../data/processed/CAFE_hist.annual.days_over_0.ehf_severity_Aus.zarr"
)
ax.set_xticks(data.lon.values, labels=data.lon.values, rotation=70)
ax.set_yticks(data.lat.values)
ax.grid()<jupyter_output><empty_output><jupyter_text>Near-surface (2m) temperature relative to AGCD<jupyter_code>plot_hindcasts(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    ["CAFE_hist", "CanESM5_hist", "EC_Earth3_hist"],
    ["AGCD", "JRA55"],
    "annual",
    "t_ref",
    "Aus_NRM",
)<jupyter_output><empty_output><jupyter_text>Anomaly correlation coefficient (Pearson)<jupyter_code>_ = plot_metrics(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    "AGCD",
    ["annual", "4-year"],
    "t_ref",
    ["rXY", "ri"],
    "Aus_NRM",
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score<jupyter_code>_ = plot_metrics(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    "AGCD",
    ["annual", "4-year"],
    "t_ref",
    ["MSSS_clim", "MSSS_pers", "MSSS_hist"],
    "Aus_NRM",
)<jupyter_output><empty_output><jupyter_text>Precipitation relative to AGCD<jupyter_code>plot_hindcasts(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    ["CAFE_hist", "CanESM5_hist", "EC_Earth3_hist"],
    ["AGCD"],
    "annual",
    "precip",
    "Aus_NRM",
)<jupyter_output><empty_output><jupyter_text>Anomaly correlation coefficient (Pearson)<jupyter_code>_ = plot_metrics(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    "AGCD",
    ["annual", "4-year"],
    "precip",
    ["rXY", "ri"],
    "Aus_NRM",
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score<jupyter_code>_ = plot_metrics(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    "AGCD",
    ["annual", "4-year"],
    "precip",
    ["MSSS_clim", "MSSS_pers", "MSSS_hist"],
    "Aus_NRM",
)<jupyter_output><empty_output><jupyter_text>Near-surface (10m) wind relative to JRA55<jupyter_code>plot_hindcasts(["CAFEf6"], ["CAFE_hist"], ["JRA55"], "annual", "V_tot", "Aus_NRM")<jupyter_output><empty_output><jupyter_text>Anomaly correlation coefficient (Pearson)<jupyter_code>_ = plot_metrics(
    ["CAFEf6"], "JRA55", ["annual", "4-year"], "V_tot", ["rXY", "ri"], "Aus_NRM"
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score<jupyter_code>_ = plot_metrics(
    ["CAFEf6"],
    "JRA55",
    ["annual", "4-year"],
    "V_tot",
    ["MSSS_clim", "MSSS_pers", "MSSS_hist"],
    "Aus_NRM",
)<jupyter_output><empty_output><jupyter_text>Extreme 2m temperature relative to AGCDHere we assess the proportion of days in forecast period with maximum temperature over 90th percentile percentile. Percentile thresholds for each day-of-year are determined for each dataset from its own climatology over the period 1991-2020. In the case of hindcast data, percentile thresholds are determined for each lead time independently. All instances of Feb 29th are removed for all datasets to simplify processing.<jupyter_code>plot_hindcasts(
    ["CAFEf6"],
    ["CAFE_hist"],
    ["AGCD"],
    "annual",
    "t_ref_max",
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Anomaly correlation coefficient (Spearman)<jupyter_code>_ = plot_metrics(
    ["CAFEf6"],
    "AGCD",
    ["annual", "4-year"],
    "t_ref_max",
    ["rXY"],
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score<jupyter_code>_ = plot_metrics(
    ["CAFEf6"],
    "AGCD",
    ["annual", "4-year"],
    "t_ref_max",
    ["MSSS_clim", "MSSS_pers", "MSSS_hist"],
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Extreme precipitation relative to AGCDHere we assess the proportion of days in forecast period with precipitation over 90th percentile percentile. Percentile thresholds for each day-of-year are determined for each dataset from its own climatology over the period 1991-2020. In the case of hindcast data, percentile thresholds are determined for each lead time independently. All instances of Feb 29th are removed for all datasets to simplify processing.<jupyter_code>plot_hindcasts(
    ["CAFEf6"],
    ["CAFE_hist"],
    ["AGCD"],
    "annual",
    "precip",
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Anomaly correlation coefficient (Spearman)<jupyter_code>_ = plot_metrics(
    ["CAFEf6"],
    "AGCD",
    ["annual", "4-year"],
    "precip",
    ["rXY"],
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score<jupyter_code>_ = plot_metrics(
    ["CAFEf6"],
    "AGCD",
    ["annual", "4-year"],
    "precip",
    ["MSSS_clim", "MSSS_pers", "MSSS_hist"],
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Extreme Forest Fire Danger Index relative to JRA55Here we assess the proportion of days in forecast period with FFDI over 90th percentile percentile. Percentile thresholds for each day-of-year are determined for each dataset from its own climatology over the period 1991-2020. In the case of hindcast data, percentile thresholds are determined for each lead time independently. All instances of Feb 29th are removed for all datasets to simplify processing.  Note that FFDI is calculated in a non-standard way here due to data availability (see Richardson et al. 2021 and Squire et al. 2021): - At the moment, the drought factor is estimated as the 20-day accumulated rainfall normalised to lie between 0 and 10, with larger values indicating less precipitation. This tends to underpredict the number of very low FFDU days, but produces similar extreme values (see `Comparison of FFDI using different drought factors`). The normalisation is carried out over all initial dates, lead times and ensemble members in the training period. This approach ensures that all days-of-year are included in the normalisation, however it does not account for lead-time dependent bias in precipitation. - We use daily mean relative humidity at 1000 hPa rather than daily max relative humidity at 2m - We use daily mean wind speed at 10m rather than daily max wind speed<jupyter_code>plot_hindcasts(
    ["CAFEf6"],
    ["CAFE_hist"],
    ["JRA55"],
    "annual",
    "ffdi",
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Anomaly correlation coefficient (Spearman)<jupyter_code>_ = plot_metrics(
    ["CAFEf6"],
    "JRA55",
    ["annual", "4-year"],
    "ffdi",
    ["rXY"],
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score<jupyter_code>_ = plot_metrics(
    ["CAFEf6"],
    "JRA55",
    ["annual", "4-year"],
    "ffdi",
    ["MSSS_clim", "MSSS_pers", "MSSS_hist"],
    "Aus_NRM",
    "days_over_p90",
)<jupyter_output><empty_output><jupyter_text>Drought relative to AGCDHere a drought event is defined as when rainfall over a 3-month period is in the lowest decile (this is the Bureau of Meteorology definition, see https://link.springer.com/article/10.1007/s00382-020-05235-0). We assess the proportion of 3-month periods in forecast period with precip within the first decile (equivalent to assessing the number of droughts in the forecast period). Note that the 3-month rainfall is not calculated in a rolling manner. Thus the max number of droughts in a year, for example, is 4. Percentile thresholds for each 3-month period are determined for each dataset from its own climatology over the period 1991-2020. In the case of hindcast data, percentile thresholds are determined for each lead time independently.<jupyter_code>plot_hindcasts(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    ["CAFE_hist", "CanESM5_hist", "EC_Earth3_hist"],
    ["AGCD"],
    "annual",
    "precip",
    "Aus_NRM",
    "3-months_under_p10",
)<jupyter_output><empty_output><jupyter_text>Anomaly correlation coefficient (Spearman)<jupyter_code>_ = plot_metrics(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    "AGCD",
    ["annual", "4-year"],
    "precip",
    ["rXY"],
    "Aus_NRM",
    "3-months_under_p10",
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score<jupyter_code>_ = plot_metrics(
    ["CAFEf6", "CanESM5", "EC_Earth3"],
    "AGCD",
    ["annual", "4-year"],
    "precip",
    ["MSSS_clim", "MSSS_pers", "MSSS_hist"],
    "Aus_NRM",
    "3-months_under_p10",
)<jupyter_output><empty_output><jupyter_text>Excess Heat Factor severity relative to AGCDHere, the excess heat factor (EHF) is calculated following Nairn and Fawcett, 2015 (https://www.mdpi.com/1660-4601/12/1/227/htm). To compute the EHF for the hindcasts, the daily temperature forecasts are first (additive) bias corrected over the verification period (1991-2020) relative to AGCD. The same is also true for the historical model data. The temperature and EHF percentiles used in the calculation of the hindcast/historical run EHF severity (T_95 and EHF_85) are calculated from AGCD over the Nairn and Fawcett periods: 1971-2000 and 1958-2011, respectively. Below we assess the proportion of days in forecast period with EHF severity over 0 (0-1: low-intensity heatwave; 1-3: severe but not extreme heatwave; > 3: extreme heatwave).We present the following on the native CAFE grid since aggregation to larger regions presents issues with interpreting the magnitude of the EHF.<jupyter_code>plot_hindcasts(
    ["CAFEf6"], ["CAFE_hist"], ["AGCD"], "annual", "ehf_severity", "Aus", "days_over_0"
)<jupyter_output><empty_output><jupyter_text>Anomaly correlation coefficient (Spearman)<jupyter_code>_ = plot_metric_maps(
    ["CAFEf6"],
    "AGCD",
    "ehf_severity",
    "rXY",
    "Aus",
    "days_over_0",
    figsize=(12, 4),
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score relative to climatology<jupyter_code>_ = plot_metric_maps(
    ["CAFEf6"],
    "AGCD",
    "ehf_severity",
    "MSSS_clim",
    "Aus",
    "days_over_0",
    vrange=(-5, 1),
    figsize=(12, 4),
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score relative to persistence<jupyter_code>_ = plot_metric_maps(
    ["CAFEf6"],
    "AGCD",
    "ehf_severity",
    "MSSS_pers",
    "Aus",
    "days_over_0",
    vrange=(-5, 1),
    figsize=(12, 4),
)<jupyter_output><empty_output><jupyter_text>Mean squared skill score relative to historical simulations<jupyter_code>EHF_fig = plot_metric_maps(
    ["CAFEf6"],
    "AGCD",
    "ehf_severity",
    "MSSS_hist",
    "Aus",
    "days_over_0",
    vrange=(-5, 1),
    figsize=(12, 4),
)<jupyter_output><empty_output><jupyter_text>Where is the MSSS positive relative to climatology, persistence and historical simulations?Lighter green indicates both are positive; darker green indicates both are positive and significant<jupyter_code>EHF_fig = plot_metric_maps(
    ["CAFEf6"],
    "AGCD",
    "ehf_severity",
    ["MSSS_clim", "MSSS_pers", "MSSS_hist"],
    "Aus",
    "days_over_0",
    figsize=(12, 4),
)<jupyter_output><empty_output><jupyter_text>Annual hindcasts for the Adelaide grid box<jupyter_code>loc = [138.75, -35.39]


def plot_loc(ax, loc):
    import cartopy.crs as ccrs

    ax.plot(
        loc[0], loc[1], "C3", marker="+", markersize=20, transform=ccrs.PlateCarree()
    )


for ax in EHF_fig.axes:
    plot_loc(ax, loc)

EHF_fig
hcst = (
    xr.open_zarr("../data/processed/CAFEf6.annual.days_over_0.ehf_severity_Aus.zarr")
    .sel(lon=loc[0], lat=loc[1], method="nearest")
    .compute()
)
hist = (
    xr.open_zarr("../data/processed/CAFE_hist.annual.days_over_0.ehf_severity_Aus.zarr")
    .sel(lon=loc[0], lat=loc[1], method="nearest")
    .compute()
)
obsv = (
    xr.open_zarr("../data/processed/AGCD.annual.days_over_0.ehf_severity_Aus.zarr")
    .sel(lon=loc[0], lat=loc[1], method="nearest")
    .compute()
)

_ = plot.hindcasts({"CAFEf6": hcst}, {"AGCD": obsv}, {"CAFE_hist": hist})<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12315959952885748, 'fraction_numerical': 0.025765606595995288, 'mean_word_length': 4.9819462791721705, 'pattern_counts': {'":': 4, '<': 121, '<?xml version=': 0, '>': 122, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '467043', 'license': 'MIT', 'n_tokens_mistral': 5120, 'n_tokens_neox': 4680, 'n_words': 1276, 'path': 'docs/results/notebooks/Australian_skill.ipynb', 'repo_name': 'dougiesquire/Squire_2022_CAFE-f6', 'size': '13584'}