Confusion about ExtremeValues adjustment

[bweeding]

Hi xclim team, I have gotten myself quite confused about the ExtremeValues function and I'm hoping for some clarity and help.

Here is an example workflow, where I'm bias adjusting daily rainfall from a CORDEX CMIP5 model over South East Asia using ERA5-LAND as the reference data. In the example I've selected a single location to make things quick and easy to plot. I've bias adjusted the daily rainfall already (monthly grouping, multiplicative quantile delta mapping, 100 quantiles) and now want to adjust the extremes.

ex_ad = xsdba.adjustment.ExtremeValues.train(ref=era5_data,hist=hist_adj_data,q_thresh=0.95,cluster_thresh='1 mm/day')

pr_hist_adj_ex = ex_ad.adjust(sim=hist_data,scen=hist_adj_data).compute()

"hist_data" is the raw model data, "hist_adj_data" has undergone bias correction.

However, when I sort each of the datasets and plot them against the ERA5 data, extreme adjust has made things worse... What am I doing wrong? Should both sim and scen in ex_ad.adjust be hist_adj_data?

Cheers
Ben

[coxipi]

Hi Ben,

I confirm I see weird things happening on my end too:

I already had the intuition that maybe we should treat this a bit differently when there is time grouping in the first adjustment (see #84 ). But even without time grouping, I get some problems and QDM+EX underperforms compared to QDM for large quantiles, and it should not.

Thanks for reporting this, I will look into this.

Best,
Éric

[bweeding]

Hi Éric, nice to know I wasn't going totally mad. I did wonder if it had something to do with the time groupings... Let me know if I can help somehow.

Cheers
Ben

[coxipi]

cluster_thresh=1 mm/d is probably too low? At least from the climates I have worked with. This is a threshold often used to determine if days are dry or wet. There are lot of days beyond this, see Vancouver from my example above

A cluster threshold of 10 mm/d has less values in comparison. Using this threshold leads to more sensible results

although Vancouver still has a weird inflexion in it. That is not the end of problems.

Choosing an appropriate threshold is not a straightforward task, see for instance:

https://revstat.ine.pt/index.php/REVSTAT/article/view/110

[bweeding]

Hmm yes I've wondered about setting a relative cluster threshold based on percentiles at each location, could do it with a for loop and dask I suppose.

I've had a play with a few different thresholds and it does tend to improve things, but as you say still isn't perfect. For example, here's an interesting one, using a threshold of 50mm:

Does a better job of correcting the very high extremes, but negatively impacts much of the data...

Also, just double checking, when I train, should the hist be raw or bias adjusted model outputs?

[coxipi]

I think the comparison of empirical quantiles is destined to be misaligned when using Extreme Values adjustment. Consider the plots above, I just focus on the left plot in my examples, I add below the result if I only do extreme adjustment (no QDM) and I still see that only-Extremes is farther away from ref than QDM in the large quantiles.

When we look at empirical quantiles, we try to compare the largest value of ref with the largest value of sim. But that is not how it works with the parametric fit. For a given value sim0, we find the related probability: P0 = GDP_sim(sim0). Then, we find what this probability would correspond to if we had the GDP params of the ref: sim0' = GDP_{ref}^{-1}(P0).

Now consider the largest quantiles of ref and sim. About ~100 of them are shown below. Comparing these quantiles, we find in my example that the CDF values of sim w.r.t. to its GDP transform are always smaller.

The values where cdf-GDP-sim = 0 means that those quantiles in sim did not qualify as extremes (either because of declustering or simply being too small, below the extreme threshold). In total, there were only 16 extremes selected in sim.
I'm getting on a tangent here, but: That kind of defeats the purpose of using Peak-Over-Threshold. In this case, we would have been better off with block-maxima, this would give us 30 years of data.

In my example, the sim is a GCM (CanESM2) and the ref is station observations (ahccd). So maybe it's not surprising the very coarse simulation yields less extremes.

Anyways, back to the main point. The last figure shows that the 16 largest precipitations in ref and sim have vastly different interpretations in term of extremes.

With the empirical quantiles way of thinking, the maximal correction would be (if we used the maximal number of quantiles):

max(ref)/max(sim) = 1.96

In practice, this is not exactly what is applied, because we simply use the correction of the last empirical quantile with a monthly grouping, but that captures the idea. Doing a bit of digging, I found that the actual adjustment factor was 1.9.

Meanwhile, with the parametric fit, the adjustment factor turned out to be much smaller

max_sim_corrected = GDP_ref_inv(GDP_sim(max_sim))
max_sim_corrected/max(sim) = 1.59

All in all, maybe a QQ-plot is not the most appropriate tool to compare extremes. We can still ask: Does the second plot indicate a problem with the data/the fitting/the parameters we have chosen? Might be the case, I'm not sure.

[bweeding]

Ah interesting, I'll have a think about the second plot!

Just to clarify, my original plot isn't of quantiles, but a scatter plot of all the data sorted...

[coxipi]

Ah yes sorry for the confusion. I plotted the same thing. I called the values in my first plot 'quantiles', but I chose as many quantiles as time steps. So it's equivalent to the ordered timeseries

[bweeding]

Haha yes!

Also I've got map_blocks working to define thresholds on each location:

def extreme_adjuster(ref_in,hist_in,hist_adj_in):

    print(ref_in.sizes)

    print(ref_in.isel(time=0).values)

    cur_thresh = str(ref_in.quantile(0.9).values)+ ' mm/day'

    print(cur_thresh)

    ex_ad = xsdba.adjustment.ExtremeValues.train(ref=ref_in,hist=hist_in,q_thresh=0.95,cluster_thresh=cur_thresh)

    out = ex_ad.adjust(sim=hist_in,scen=hist_adj_in).transpose('time','lat','lon')

    print(out.sizes)

    return out

then

result = xr.map_blocks(
    extreme_adjuster,
    era5_data,
    args=[
    hist_data,
    hist_adj_data],
    template = hist_adj_data
).compute()

Output shows different cur_thresh values for each spatial location:

Frozen({'time': 13140, 'lat': 1, 'lon': 1})
[[5.84357977]]
26.432598495483404 mm/day
Frozen({'time': 13140, 'lat': 1, 'lon': 1})
[[5.93152618]]
24.175616264343265 mm/day
Frozen({'time': 13140, 'lat': 1, 'lon': 1})
[[8.88751411]]
22.391751480102542 mm/day
Frozen({'time': 13140, 'lat': 1, 'lon': 1})
[[7.93474913]]
24.410479164123537 mm/day
Frozen({'time': 13140, 'lat': 1, 'lon': 1})
Frozen({'time': 13140, 'lat': 1, 'lon': 1})
Frozen({'time': 13140, 'lat': 1, 'lon': 1})
Frozen({'time': 13140, 'lat': 1, 'lon': 1})

[coxipi]

thresholds on each location

Very nice! Hmm we do have a way to implement this, let me see.

EDIT: see #179

[bweeding]

Ah cool! So you will be able to generate a data array that specifies the cluster_thresh at each location?

[coxipi]

Exact, you can now pass cluster_thresh as a DataArray (I just merged #179 and we will also have a release of xsdba in the coming weeks)!

[bweeding]

Great, that should make investigating the issues much easier!

[bweeding]

Is the release imminent or should I try working with a zip of the current package?

[Zeitsperre]

Hi @bweeding We should be releasing a new version in the coming days!

[Zeitsperre]

This has been released. Should be on conda-forge by tomorrow!