pangeo-data / xESMF

Universal Regridder for Geospatial Data
http://xesmf.readthedocs.io/
MIT License
182 stars 32 forks source link

Dimension order shouldn't matter for periodicity #363

Open mberdahl-uw opened 1 month ago

mberdahl-uw commented 1 month ago

Hi,

I'm trying to regrid some cmip6 models with xesmf regridder. The operation works, but upon plotting some (not all) of the regridded models have a seam in the output at the 0 dateline. I am using periodic=True. Appreciate any help!

Here is a small bit of code that illustrates the issue with one model: Data.zip

import xarray as x
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import numpy as np
import xesmf as xe

original_dataset = xr.open_dataset("ds_input.nc")
ds_target = xr.open_dataset("ds_target.zarr")

# Define the names of the datasets
dataset_names = ["original_dataset"]

# Create an empty dictionary to store the regridded datasets
regridded_datasets = {}

# Loop through each dataset name
for dataset_name in dataset_names:

    # Define the original dataset
    original_dataset = globals()[dataset_name]

    # Define the regridder object
    regridder = xe.Regridder(
        original_dataset, ds_target, "bilinear", periodic=True, ignore_degenerate=True,
    )  # this takes some time to calculate a weight matrix for the regridding

    # Apply regridding and store the regridded dataset in the dictionary
    regridded_datasets[dataset_name] = regridder(original_dataset.squeeze())  # this is a bit of a lazy operation

original_regridded = regridded_datasets["original_dataset"]

# for testing, just compare the original with the regridded dataset to see about the seam

fig, (axs) = plt.subplots(
    # ncols=2, nrows=11, figsize=[10, 35], subplot_kw={"projection": ccrs.Robinson()}
    ncols=2, nrows=1, figsize=[10, 3], subplot_kw={"projection": ccrs.PlateCarree()}

)

vupper = 293.15
vlower = 253.15

# plot the model data
original_dataset.tas.plot(
    ax=axs[0],
    x="x",
    y="y",
    transform=ccrs.PlateCarree(),
    cmap="RdBu_r",
    vmax = vupper,
    vmin = vlower,
    robust=True,
)

original_regridded.tas.plot(
    ax=axs[1],
    x="x",
    y="y",
    transform=ccrs.PlateCarree(),
    cmap="RdBu_r",
    vmax = vupper,
    vmin = vlower,
    robust=True,
)

axs[0].coastlines()
axs[0].set_title("Historical ACCESS_CM2")

axs[1].coastlines()
axs[1].set_title("Historical ACCESS_CM2 - regridded")

Screenshot 2024-05-15 at 10 55 46 AM (2)

image

aulemahal commented 1 month ago

Hi @mberdahl-uw, thanks for the very complete issue.

The bug is quite subtle and needs some fix as fast as possible... It has to do with the dimension order of your data. Your ds_input has a lon coordinate of dimensions x, y, while ESMF expects y, x. xESMF should have taken care of this, dimension order is not supposed to be a problem in the xarray-world. But it has failed to do so and for some reason this has made the periodic feature fail, while the rest of the regridding seems to work anyway.

On my side, I was able to produce the correct figures by adding original_dataset = original_dataset.transpose('y', 'x') before initialising the Regridder.

Please confirm if this fixes your problem, we'll work on a way to fix this.

mberdahl-uw commented 1 month ago

Hi @aulemahal , thanks so much for looking into this. Your fix works for me too, so that is great. I'll look out for a fix on your end when you get around to it. In the meantime I'll do this piecemeal approach for the models that require this transpose.