Empirical Quantile Mapping - Implementation

How EQM is computed in src/bias_correction.py and src/distribution_fitting.py.

For the theory, see Methodology > EQM. This page is the algorithm view.

Entry point

EQM runs after LS, inside the same lseqm() function in src/bias_correction.py. The per-cell quantile mapping is delegated to gamma_quantile_mapping_precomputed() in src/distribution_fitting.py.

Per-cell loop

EQM is a per-pixel operation. The pipeline iterates (i, j) over the AOI grid:

for i in range(n_lat):
    for j in range(n_lon):
        imerg_ts = ls_corrected_precip.values[:, i, j]
        if np.all(np.isnan(imerg_ts)):
            continue   # ocean - skip

        if _use_native_cpc:
            # Use pre-computed CPC parameters (smooth path)
            eqm_data[:, i, j] = gamma_quantile_mapping_precomputed(
                imerg_ts,
                cpc_gamma_shape   = interp_cpc_params['gamma_shape'][i, j],
                cpc_gamma_scale   = interp_cpc_params['gamma_scale'][i, j],
                cpc_gpd_threshold = interp_cpc_params['gpd_threshold'][i, j],
                ... # all CPC params already interpolated to (i, j)
            )
        else:
            # Per-pixel fitting (no native CPC path)
            cpc_ts = cpc_dekad_data.values[:, i, j]
            eqm_data[:, i, j] = gamma_quantile_mapping(imerg_ts, cpc_ts)

The three-distribution model

gamma_quantile_mapping_precomputed actually fits three distributions for CPC:

  1. Dry-day mass at precip < WET_DAY_THRESHOLD (default 1 mm/day). Probability p_dry_cpc.
  2. Gamma body for the wet-day distribution below the GPD threshold.
  3. GPD tail above the 80th percentile of wet days (see GPD Tail Implementation).

The IMERG side gets its own Gamma fit on the wet-day sample for the dekad. The unconditional CDF is then:

F_uncond(x) = p_dry + (1 - p_dry) * F_wet(x)

and the same form is used to invert through CPC’s three-distribution model.

This is Cannon-style dry-day handling combined with a parametric tail correction. The full formula and derivation is in the docstring of gamma_quantile_mapping_precomputed.

Dry-day handling (correct-zero)

Two distinct paths produce a zero output, both inside gamma_quantile_mapping / gamma_quantile_mapping_precomputed:

  1. Dry input stays dry. The output array is zero-initialised (out = np.zeros_like(...)), and only above-threshold IMERG days are remapped (out[is_wet_imerg] = out_wet). A day the satellite already reports below WET_DAY_THRESHOLD is never touched and emerges as exactly 0.0.
  2. Over-detected drizzle is killed. The Cannon 2015 §3.2 step builds the unconditional CDF position of each wet IMERG value, y_uncond = p_dry_imerg + (1 - p_dry_imerg) * F_wet(x), and sets to zero every value whose position falls below the reference dry-day fraction: kill = y_uncond < p_dry_cpc. This is what matches the corrected wet-day frequency to CPC.
out = np.zeros_like(imerg_valid)            # dry days and killed drizzle -> 0
y_uncond = p_dry_imerg + (1 - p_dry_imerg) * gamma.cdf(imerg_wet, shape1, scale=scale1)
kill = y_uncond < p_dry_cpc                 # over-detected drizzle
...
out_wet = np.where(kill, 0.0, out_wet)
out[is_wet_imerg] = out_wet                 # only wet inputs are remapped

p_dry_cpc (the reference dry-day fraction) is computed in fit_cpc_parameters_on_native_grid as 1.0 - len(wet)/len(valid). LS is multiplicative, so it preserves exact zeros on its own (0 * factor = 0); the reclassification of over-detected drizzle is entirely the EQM stage, and the CNN refinement leaves dry days untouched. See Methodology > Wet / dry handling for the verified dry-day agreement numbers.

Where CPC parameters come from

Two code paths, picked by whether cpc_native_ds is passed:

Path CPC params source
Smooth (BCSD) fit_cpc_parameters_on_native_grid(cpc_native_dekad) -> dict of DataArrays at 0.5 deg, then interpolate_cpc_params_to_imerg_grid(...) -> dict of DataArrays at 0.1 deg. Bilinear + reindex(nearest) fillna.
Per-pixel Fit on the fly inside gamma_quantile_mapping() for each (i, j) cell using cpc_ts. Slower; visible block-boundary artefact at the original CPC 0.5 deg resolution.

The notebooks default to the smooth path. The smooth (BCSD) path is what the figure in Methodology > BCSD illustrates: parameters fit blocky at 0.5 deg, bilinearly disaggregated to smooth 0.1 deg, then per-pixel quantile mapping.

Cost

For Indonesia (171 x 461 cells, ~25 years x 10 days per dekad = ~250 daily samples per cell):

  • Smooth path: dominated by the 79k cells x ~10 ms quantile mapping = ~13 min per dekad. Times 36 dekads = ~8 hours of nb02 compute on CPU.
  • Per-pixel path: adds ~2-3 ms per cell for the Gamma fit; ~25 min per dekad.

Bali (9 x 14 = 126 cells) is sub-second per dekad either way.

Failure modes and skips

  • All-NaN IMERG cells (ocean) are skipped (continue at line 211 in bias_correction.py).
  • Cells where any CPC parameter is NaN return all-NaN (line 528-530 in distribution_fitting.py).
  • Cells with fewer than 5 wet days return all-zero (line 549-553). This is conservative - the alternative is to propagate IMERG raw, which may be wetter than the climatology supports.
Back to top