Fix NaN propagation in ISIMIP v3.0 non-parametric quantile mapping (_isimip_quantile_map_x_on_y_non_parametically)

[AntonioCC808]

Description

Currently, the ISIMIP v3.0 non-parametric quantile mapping helper _isimip_quantile_map_x_on_y_non_parametically (located in ibicus/utils/_math_utils.py) is vulnerable to missing values (NaN). A single NaN in a grid point's time-series window causes the entire window to be filled with NaNs in the debiased output.

The Issues

1. Default nan_policy in rankdata: Without specifying nan_policy in scipy.stats.rankdata(x), it defaults to propagating NaNs. If x has a single NaN, all ranks become NaN, resulting in p_x and the output z being entirely NaNs.
2. Denominator Size: Using x.size as the denominator when x contains NaNs shifts the empirical probabilities down.
3. Reference array y NaNs: If the reference/observations array y contains NaNs, np.sort(y) places them at the end. Since p_y = np.linspace(0.0, 1.0, y.size) is mapped to the sorted array, np.interp will return NaN for any high value in x (specifically whenever p_x > valid_fraction), causing artificial loss of valid data in the upper tail of the distribution.

Proposed Solution

Modify _isimip_quantile_map_x_on_y_non_parametically to:

1. Ignore NaNs when calculating ranks in x using nan_policy='omit'.
2. Use the count of non-NaN elements in x (valid_count) as the denominator.
3. Filter out NaNs from the reference array y before sorting and mapping (y_valid = y[~np.isnan(y)]). This prevents np.interp from propagating NaNs from the upper tail of y into valid outputs of z.
4. Include a safety check for windows where all values are NaN (valid_count == 0 or y_valid.size == 0) to prevent division by zero or errors in np.interp.

Proposed Code Changes

def _isimip_quantile_map_x_on_y_non_parametically(
    x: np.ndarray, y: np.ndarray
) -> np.ndarray:
    # 1. Filter out NaNs from reference observations to prevent interpolating against NaNs
    y_valid = y[~np.isnan(y)]
    
    # 2. Compute ranks omitting NaNs, keeping NaNs in their original positions
    ranks = scipy.stats.rankdata(x, nan_policy='omit')
    
    # 3. Denominator should be the count of valid elements in x, not x.size
    valid_count = np.count_nonzero(~np.isnan(x))
    
    # 4. Guard check for empty/all-NaN inputs to avoid division by zero or empty interpolation
    if valid_count == 0 or y_valid.size == 0:
        return np.full_like(x, np.nan, dtype=y.dtype)
        
    # 5. Compute empirical ranks probability
    p_x = (ranks - 1.0) / valid_count
    
    p_y = np.linspace(0.0, 1.0, y_valid.size, dtype=y.dtype)
    z = np.interp(p_x, p_y, np.sort(y_valid)).astype(y.dtype)
    return z

[jakobwes]

Hi, thank you for creating an issue and proposing a fix. This is extremely helpful!

We had a brief look at the issue. Whilst, generally we believe that better NaN handling would be very valuable within ibicus, I am a bit wary of doing this for individual methods in an inconsistent way. For example, even with the above fix, ISIMIP would still fail in the presence of NaNs when parametric quantile mapping is used, as the scipy distribution fits doesn't handle them, I believe. If failsafe = True is activated this would still lead to the entire window being filled with NaNs. A more principled way would probably be to add NaN infilling or an ignore_nan-argument to the debiaser apply function which could filter out NaNs before any bias adjustment mapping. However, we will have to discuss this internally as these are also things that could be done on the user side and there are limits to how many modifications and checks on user data we would like to make within ibicus: the package is meant to operate purely on numerical data some amount of trust in the user formatted data is necessary.

We will discuss within the maintainer team. Please do let us know however if NaN handling is crucial for your use case and whether filtering them out before giving the data to ibicus would cause issues.

Best,
Jakob