Clarifying instructions in "First example : pr and frequency adaptation" example

[bweeding]

Generic Issue

Posting as I'm confused about the example code provided for adjusting precipitation data including frequency adjustment.

In the code below from https://xclim.readthedocs.io/en/stable/notebooks/sdba.html#First-example-:-pr-and-frequency-adaptation, the training appears to take place without reference to pr_hist:

# 2nd try with adapt_freq
sim_ad, pth, dP0 = sdba.processing.adapt_freq(
    pr_ref, pr_sim, thresh="0.05 mm d-1", group="time"
)
QM_ad = sdba.EmpiricalQuantileMapping.train(
    pr_ref, sim_ad, nquantiles=15, kind="*", group="time"
)
scen_ad = QM_ad.adjust(pr_sim)

Is that correct? Or should pr_hist be involved in generating QM_ad, which is then applied to sim_ad?

Cheers
Ben

Code of Conduct

• I agree to follow this project's Code of Conduct

[aulemahal]

Oh, hey you are right, this is an error! Indeed, pr_hist should be the second input of adapt_freq and hist_ad, its output, be the input of the train step.

However, pr_sim is correctly the input of the adjust step here. The primary goal of adapt_freq is to avoid problems in the adjustment factors computed by the train. Once we have the bug-free adjustment factors, there's no reason to use the inflated version of the precipitation to create scen.

At least, that's what I understand from the reference article. Their adjustment method is not exactly xclim's EQM, but the idea is the same.
Themeßl et al. - 2012 - Empirical-statistical downscaling and error correc.pdf

[bweeding]

Ah good, that's what I've done in my code:

hist_fad, pth, dP0 = sdba.processing.adapt_freq(
    agcd.pr.chunk(chunks={"time":10950,"lon": 25, "lat": 25}), cnrm_hist.pr.chunk(chunks={"time":10950,"lon": 25, "lat": 25}), thresh="0.05 mm d-1", group="time"
)

proj_fad, pth, dP0 = sdba.processing.adapt_freq(
    agcd.pr.chunk(chunks={"time":10950,"lon": 25, "lat": 25}), cnrm_proj85.pr.chunk(chunks={"time":34310,"lon": 25, "lat": 25}), thresh="0.05 mm d-1", group="time"
)

cur_trainer = sdba.QuantileDeltaMapping.train(
    agcd.pr.chunk(chunks={"time":10950,"lon": 25, "lat": 25}), hist_fad.chunk(chunks={"time":10950,"lon": 25, "lat": 25}), nquantiles=15, kind="*", group="time"
)

hist_adj_out = cur_trainer.adjust(cnrm_hist.pr.chunk(chunks={"time":10950,"lon": 25, "lat": 25}), interp="linear")

proj_adj_out = cur_trainer.adjust(cnrm_proj85.pr.chunk(chunks={"time":34310,"lon": 25, "lat": 25}), interp="linear")