Properly handling of non permanent masks

[stefraynaud]

It happens that gridded data have missing values at locations that are not permanently masked.
One way to handle this is to regrid the mask it self, once converted to float. Then, the final mask should depend of the regridded mask values and on the regridding method.

At least an option should be added to handle this situation.

[stefraynaud]

Here is an example of how an adaptative masking is applicable to data with non permanent missing values, when using a conservative method: this helps not masking cell that are only partially masked.

Classic case:

Adaptative masking:

Here is the code:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Oct 21 09:24:15 2020 by sraynaud
"""
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import xesmf as xe


def double_plot(da0, da1, suptitle, figfile):
    fig, (ax0, ax1) = plt.subplots(ncols=2, sharex=True, sharey=True,
                                   subplot_kw={'aspect': 1}, figsize=(8, 4))
    da0.plot(ax=ax0, add_colorbar=False, vmin=0, vmax=1)
    ax0.set_title('Input')
    ax0.grid()
    da1.plot(ax=ax1, add_colorbar=False, vmin=0, vmax=1)
    ax1.set_title('Output')
    ax1.grid()
    plt.axis([0, 3, 0, 3])
    plt.suptitle(suptitle, weight="bold")
    plt.savefig(figfile)
    plt.show()
    plt.close()


# %% Input dataset
tempi = np.arange(9).reshape(3, 3)*1.*0+1
tempi[1, 1] = np.nan
dsi = xr.Dataset(
    {'temp': (['lat', 'lon'], tempi),
     'lon_b': ('nxb', [0., 1., 2., 3.]),
     'lat_b': ('nyb', [0., 1., 2., 3.])},
    coords={'lat': [.5, 1.5, 2.5], 'lon': [.5, 1.5, 2.5]})

# %% Output grid
grido = xr.Dataset(
    {'lon_b': ('nxb', [.25, 0.75, 1.25, 1.75, 2.25, 2.75]),
     'lat_b': ('nyb', [.25, 0.75, 1.25, 1.75, 2.25, 2.75])},
    coords={'lat': [.5, 1., 1.5, 2., 2.5], 'lon': [.5, 1., 1.5, 2., 2.5]})

# %% Regridder
rg = xe.Regridder(dsi, grido, method="conservative")

# %% Regridding
tempo = rg(dsi['temp'])

# %% Plots
double_plot(dsi['temp'], tempo, "Without adaptative masking (classic)",
            __file__[:-2]+'default.png')

# %% Adaptative masking
validi = dsi['temp'].notnull().astype('d')
valido = rg(validi)
tempi0 = dsi['temp'].fillna(0)
tempo0 = rg(tempi0)
tempom = xr.where(valido != 0, tempo0 / valido, np.nan)
double_plot(dsi['temp'], tempom, "With adaptative masking",
            __file__[:-2]+'adaptative_masking.png')

[sol1105]

In the following jupyter notebook I took a closer look into the adaptive masking technique and in regridding with the bilinear and the patch method. I found there to be a weird behaviour at the southern domain boundary for those two methods, as they introduce missing values. I would be interested in your opinion on this behaviour and in general if it makes sense to apply the adaptive masking technique for other methods than the conservative method.
https://github.com/roocs/regrid-prototype/blob/main/docs/notebooks/xESMF_nonpermanent_masks.ipynb

In another notebook I then applied the adaptive masking technique on low and high resolution MPI-ESM1-2 MPIOM output as published to CMIP6, which might be also interesting to look at:
https://github.com/roocs/regrid-prototype/blob/main/docs/notebooks/xESMF_test_Adaptive_Masking-MPIOM.ipynb

[sol1105]

Salut Stephane,

I updated the second notebook (the one with the application of the adaptive masking on output of an ocean model):
https://github.com/roocs/regrid-prototype/blob/main/docs/notebooks/xESMF_test_Adaptive_Masking-MPIOM.ipynb

For the adaptive masking I suggest to define a minimum contribution threshold as follows:

An addition to the adaptive masking / renormalization method would be to apply a minimum threshold of contributing source cell area to the target grid cell. In this way one could allow the renormalization only if at least for example 80% of the source cells area that contributes to a target cell is not masked. Else the result would be a missing value.

def adaptive_masking(ds_in, regridder, min_norm_contribution=1):
    """Performs regridding incl. renormalization for conservative weights"""    
    validi = ds_in.notnull().astype('d')
    valido = regridder(validi)
    tempi0 = ds_in.fillna(0)
    tempo0 = regridder(tempi0)
    # min_norm_contribution factor could prevent values for cells that should be masked.
    # It prevents the renormalization for cells that get less than min_norm_contribution 
    #  from source cells. If the factor==0.66 it means that at most one third of the source cells' area
    #  contributing to the target cell is masked. This factor has however to be tweaked manually for each 
    #  pair of source and destination grid.
    if min_norm_contribution<1:
        valido = xr.where(valido < min_norm_contribution, np.nan, valido)
    ds_out = xr.where(valido != 0, tempo0 / valido, np.nan)
    return ds_out

I would like your feedback on that and the few questions posed in the notebook.

Thank you and kind regards
Martin

[stefraynaud]

Ok, I'll have a look at it soon.

[stefraynaud]

@sol1105 That's a really nice work.
It is a good point that we find similar results between without_adaptative_mask+add_matrix_nans and with_adaptative_masking.
It seems that adapative masking is mandatory for the patch method.
The min_norm_contribution is very useful addition, especially for land-sea masking.

Here is a possible strategy:

• Apply the add_matrix_nans in all cases.
• Check if the user want the adaptative_masking (or let it True my default? I'm more in favour of this solution) or check if the patch method is used.
• If adaptative masking is desired or needed, then check if the mask is permanent and compare it to the provided input mask if any, then decide if adaptative masking is necessary.
• Apply it.

[sol1105]

I applied the adaptive masking for the patch and bilinear methods on the MPI-ESM1-2 MPIOM model output and updated the notebook accordingly (the new bits are at the bottom): https://github.com/roocs/regrid-prototype/blob/main/docs/notebooks/xESMF_test_Adaptive_Masking-MPIOM.ipynb

This might help deciding for which methods (besides conservative) the adaptive masking should be enabled by default. For the patch method it looks advantageous to use the technique, but I would not use it for the bilinear method (at least not without any additional settings like a target mask or a mask-threshold, which I did not test however).

[stefraynaud]

Maybe the default masking threshold should be set to 0.5. This would make the destination mask closer in shape to the source mask for all methods.