xCDAT / xcdat

An extension of xarray for climate data analysis on structured grids.
https://xcdat.readthedocs.io/en/latest/
Apache License 2.0
101 stars 11 forks source link

[Bug]: NaN handling in regrid2 #580

Closed lee1043 closed 3 months ago

lee1043 commented 7 months ago

What happened?

When NaN value exists in the original field, interpolation using regrid2 replaces NaN values to zero, then interpolates, which introduces error in the resulted field. unmapped_to_nan=True that was discussed in #528 for xesmf is not available for regrid2, wondering if that could be added.

What did you expect to happen? Are there are possible answers you came across?

output

Original field has NaN values over the ocean, and when interpolate, it should look like one in the middle (example used xesmf). However, regrid2 injects zeros to NaN values before interpolation, then interpolates, resulting zero over land and gradation values along coastlines.

Minimal example code to reproduce the above plot is attached below.

Minimal Complete Verifiable Example (MVCE)

import matplotlib.pyplot as plt
import xcdat as xc

# Generate a sample data field that has value only over ocean
input_data = "ts_Amon_ACCESS1-0_historical_r1i1p1_185001-200512.nc"
input_data2 = "sftlf_fx_ACCESS1-0_amip_r0i0p0.nc"

ds = xc.open_dataset(input_data)
ds2 = xc.open_dataset(input_data2)
sftlf = ds2["sftlf"]
ds["sst"] = ds["ts"].where(sftlf==0)

# Make a target grid to interpolate to
target_grid = xc.create_uniform_grid(-90, 90, 5, 2.5, 357.5, 5)

# regrid
ds_regridded_xesmf = ds.regridder.horizontal("sst", target_grid, tool="xesmf", method="bilinear")
ds_regridded_regrid2 = ds.regridder.horizontal("sst", target_grid, tool="regrid2") # unmapped_to_nan=True is not available for regrid2

fig, axs = plt.subplots(1, 3, figsize=(15, 4))

ds["sst"].isel(time=0).plot(ax=axs[0])
ds_regridded_xesmf["sst"].isel(time=0).plot(ax=axs[1])
ds_regridded_regrid2["sst"].isel(time=0).plot(ax=axs[2])

axs[0].set_title("Original")
axs[1].set_title("Regridded: xesmf")
axs[2].set_title("Regridded: regrid2")

fig.tight_layout()

Relevant log output

No response

Anything else we need to know?

Sample input files used above:

Environment

xcdat 0.6.1

xr.show_versions() INSTALLED VERSIONS ------------------ commit: None python: 3.10.10 | packaged by conda-forge | (main, Mar 24 2023, 20:17:34) [Clang 14.0.6 ] python-bits: 64 OS: Darwin OS-release: 22.6.0 machine: x86_64 processor: i386 byteorder: little LC_ALL: None LANG: None LOCALE: (None, 'UTF-8') libhdf5: 1.14.2 libnetcdf: 4.9.2 xarray: 2023.11.0 pandas: 2.1.3 numpy: 1.23.5 scipy: 1.11.4 netCDF4: 1.6.5 pydap: None h5netcdf: None h5py: None Nio: None zarr: None cftime: 1.6.3 nc_time_axis: 1.4.1 iris: None bottleneck: None dask: 2023.12.0 distributed: 2023.12.0 matplotlib: 3.7.1 cartopy: 0.22.0 seaborn: 0.12.2 numbagg: None fsspec: 2023.12.0 cupy: None pint: None sparse: 0.14.0 flox: None numpy_groupies: None setuptools: 67.7.2 pip: 23.1.2 conda: None pytest: None mypy: None IPython: 8.18.1 sphinx: None
lee1043 commented 7 months ago

Hi @jasonb5, is there any chance this could be fixed in #533?

jasonb5 commented 6 months ago

@lee1043 This will be resolved in another PR after #533, it needs to be fixed for all the regridding tools.

lee1043 commented 3 months ago

Close this issue as it was resolved by #533.

Result of the minimal example code above after updating the xcdat to the latest main branch. output

tomvothecoder commented 3 months ago

Close this issue as it was resolved by #533.

Result of the minimal example code above after updating the xcdat to the latest main branch. output

Thank you for confirming!