<jupyter_start><jupyter_text>The frequency of severe EHF (> 1) events is overestimated in the CAFEf6 forecasts when the EHF is calculated from bias corrected daily temperature. This notebook tries to understand why.Note that using mutliplicative bias correction of T, rather than additive, makes negligible difference to the following<jupyter_code>import xarray as xr

from src import utils

import matplotlib.pyplot as plt
%load_ext autoreload
%autoreload 2
%load_ext lab_black<jupyter_output><empty_output><jupyter_text>Cluster<jupyter_code>from distributed import Client
from dask_jobqueue import PBSCluster

walltime = "02:00:00"
cores = 48
memory = "192GB"
cluster = PBSCluster(
    processes=1,
    walltime=str(walltime),
    cores=cores,
    memory=str(memory),
    job_extra=[
        "-q express",
        "-l ncpus=" + str(cores),
        "-l mem=" + str(memory),
        "-P xv83",
        "-l jobfs=100GB",
        "-l storage=gdata/xv83+gdata/oi10+gdata/ua8",
    ],
    local_directory="$PBS_JOBFS",
    header_skip=["select"],
)

cluster.scale(jobs=1)
client = Client(cluster)
print(f"Dask daskboard link: {client.dashboard_link}")<jupyter_output>/g/data/xv83/users/ds0092/software/miniconda3/envs/forecast_analysis/lib/python3.9/site-packages/dask_jobqueue/core.py:20: FutureWarning: tmpfile is deprecated and will be removed in a future release. Please use dask.utils.tmpfile instead.
  from distributed.utils import tmpfile<jupyter_text>First, let's check our bias correction<jupyter_code>f6 = xr.open_zarr("../../data/processed/CAFEf6.daily.abc_1991-2020.t_ref_Aus.zarr")
agcd = xr.open_zarr("../../data/processed/AGCD.daily.full.t_ref_Aus.zarr")
period = ["1991-01-01", "2020-12-31"]

agcd_clim = (
    utils.keep_period(agcd, period).groupby("time.dayofyear").mean("time").compute()
)
f6_clim = utils.keep_period(f6, period).mean(["init", "member"]).compute()
lead = 500

f6_clim_at_lead = f6_clim.isel(lead=lead)
agcd_clim_at_lead = agcd_clim.sel(dayofyear=((lead + 305) % 365))
(f6_clim_at_lead - agcd_clim_at_lead)["t_ref"].plot()<jupyter_output><empty_output><jupyter_text>How do individual timeseries look?<jupyter_code>init = 20
sel = dict(lat=10, lon=10)

f6_ts = f6["t_ref"].isel({**sel, "init": init}).swap_dims({"lead": "time"}).compute()
agcd_ts = agcd["t_ref"].isel(sel).sel(time=f6_ts.time).compute()
fig = plt.figure(figsize=(15, 5))

plt.fill_between(
    f6_ts.time.values,
    f6_ts.min("member"),
    f6_ts.max("member"),
    color="C0",
    alpha=0.5,
    edgecolor=None,
)
plt.fill_between(
    f6_ts.time.values,
    f6_ts.quantile(0.1, "member"),
    f6_ts.quantile(0.9, "member"),
    color="C2",
    alpha=0.7,
    edgecolor=None,
)
agcd_ts.plot(color="k")
f6_ts.isel(member=0).plot(color="C1")
f6_ts.mean("member").plot(color="C3")

plt.title("Additive bias correction")<jupyter_output><empty_output><jupyter_text>What about for EHF?<jupyter_code>f6["ehf"] = utils.calculate_EHF(
    f6,
    T_p95_file="data/processed/AGCD.daily.p95_1971-2000.t_ref_Aus.zarr",
    rolling_dim="lead",
)["ehf"]
agcd["ehf"] = xr.open_zarr("../../data/processed/AGCD.daily.full.ehf_Aus.zarr")["ehf"]
init = 20
sel = dict(lat=10, lon=10)

f6_ts = f6["ehf"].isel({**sel, "init": init}).swap_dims({"lead": "time"}).compute()
agcd_ts = agcd["ehf"].isel(sel).sel(time=f6_ts.time).compute()
fig = plt.figure(figsize=(15, 5))

plt.fill_between(
    f6_ts.time.values,
    f6_ts.min("member"),
    f6_ts.max("member"),
    color="C0",
    alpha=0.5,
    edgecolor=None,
)
plt.fill_between(
    f6_ts.time.values,
    f6_ts.quantile(0.1, "member"),
    f6_ts.quantile(0.9, "member"),
    color="C2",
    alpha=0.7,
    edgecolor=None,
)
agcd_ts.plot(color="k")
f6_ts.isel(member=0).plot(color="C1")
f6_ts.mean("member").plot(color="C3")<jupyter_output>/g/data/xv83/users/ds0092/software/miniconda3/envs/forecast_analysis/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1395: RuntimeWarning: All-NaN slice encountered
  result = np.apply_along_axis(_nanquantile_1d, axis, a, q,<jupyter_text>And EHF severity?<jupyter_code>f6["ehf_severity"] = xr.open_zarr(
    "../../data/processed/CAFEf6.daily.full.ehf_severity_Aus.zarr"
)["ehf_severity"]
agcd["ehf_severity"] = xr.open_zarr(
    "../../data/processed/AGCD.daily.full.ehf_severity_Aus.zarr"
)["ehf_severity"]
init = 20
sel = dict(lat=10, lon=10)

f6_ts = (
    f6["ehf_severity"].isel({**sel, "init": init}).swap_dims({"lead": "time"}).compute()
)
agcd_ts = agcd["ehf_severity"].isel(sel).sel(time=f6_ts.time).compute()
fig = plt.figure(figsize=(15, 5))

plt.fill_between(
    f6_ts.time.values,
    f6_ts.min("member"),
    f6_ts.max("member"),
    color="C0",
    alpha=0.5,
    edgecolor=None,
)
plt.fill_between(
    f6_ts.time.values,
    f6_ts.quantile(0.1, "member"),
    f6_ts.quantile(0.9, "member"),
    color="C2",
    alpha=0.7,
    edgecolor=None,
)
agcd_ts.plot(color="k")
f6_ts.isel(member=0).plot(color="C1")
f6_ts.mean("member").plot(color="C3")<jupyter_output>/g/data/xv83/users/ds0092/software/miniconda3/envs/forecast_analysis/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1395: RuntimeWarning: All-NaN slice encountered
  result = np.apply_along_axis(_nanquantile_1d, axis, a, q,<jupyter_text>This might look somewhat reasonable, but the actually number of heatwave days predicted for each member (e.g. orange line) is substantially higer (~ 2x) what is observed:<jupyter_code>(agcd_ts > 0).sum("time").values
(f6_ts > 0).sum("time").values<jupyter_output><empty_output><jupyter_text>Thus the proportion per year is too high<jupyter_code>f6_p = xr.open_zarr(
    "../../data/processed/CAFEf6.annual.days_over_0.ehf_severity_Aus.zarr"
)
agcd_p = xr.open_zarr(
    "../../data/processed/AGCD.annual.days_over_0.ehf_severity_Aus.zarr"
)
lead = 1
sel = dict(lat=10, lon=10)

f6_ts = utils.keep_period(
    f6_p["ehf_severity"]
    .isel({**sel, "lead": lead})
    .swap_dims({"init": "time"})
    .compute(),
    period,
)
agcd_ts = agcd_p["ehf_severity"].isel(sel).sel(time=f6_ts.time).compute()
fig = plt.figure(figsize=(15, 5))

plt.fill_between(
    f6_ts.time.values,
    f6_ts.min("member"),
    f6_ts.max("member"),
    color="C0",
    alpha=0.4,
    edgecolor=None,
)
plt.fill_between(
    f6_ts.time.values,
    f6_ts.quantile(0.1, "member"),
    f6_ts.quantile(0.9, "member"),
    color="C2",
    alpha=0.6,
    edgecolor=None,
)
agcd_ts.plot(color="k")
f6_ts.isel(member=10).plot(color="C1")
f6_ts.mean("member").plot(color="C3")
cluster.close()
client.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1743810548977395, 'fraction_numerical': 0.040135322159003534, 'mean_word_length': 4.885972850678733, 'pattern_counts': {'":': 8, '<': 30, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '204081', 'license': 'MIT', 'n_tokens_mistral': 2840, 'n_tokens_neox': 2676, 'n_words': 379, 'path': 'notebooks/exploratory/explore_severe_EHF.ipynb', 'repo_name': 'dougiesquire/Squire_2022_CAFE-f6', 'size': '6503'}