Is xesmf.util.grid_global() leaving out data at poles?

[trimpinm]

I am attempting to use the grid_global function to define an output grid for global regridding. However, upon successful regridding it appears as though there is no data at the -90 and 90 degree latitude lines. I am not sure why. I have tried a variety of different resolutions, but none of them are able to preserve the data values. I am using the attached dataset, and I have included screenshots pre and post-regrid. grid_global supp.zip

[sol1105]

Hi @mtrimp2,
xesmf.util.grid_global() does not leave out data at the poles, however, the poles are only included in the bounds and the north/southmost cell centers are therefore not 90.0 and -90.0. For example, creating a global grid with 1 degree resolution will yield bounds of [-90, -89] and cell centers at -89.5 for the soutmost row of cells.

I tried it out myself and the resulting plot included the polar rows, when using the CF-conformal representation of the target grid. It seems that the plot program behaves differently for 1D and 2D coordinate variables:

import numpy as np
import xarray as xr
import xesmf as xe
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

ds = xr.open_dataset("grass_MERRA2_400.tavg1_2d_slv_Nx.20200101_variable.nc4")

# Create weights for global target grid (in CF and non-CF representation)
regridder = xe.Regridder(ds, xe.util.grid_global(1., 1.), "conservative")
regridder_cf = xe.Regridder(ds, xe.util.grid_global(1., 1., cf=True), "conservative")

# Regrid
ds_out = regridder(ds.H500)
ds_out_cf = regridder_cf(ds.H500)

# Assert both results are equal
np.testing.assert_array_equal(ds_out.values, ds_out_cf.values)

# Plot
fig, axes = plt.subplots(nrows=3, figsize=(18,12), subplot_kw={'projection': ccrs.PlateCarree()});
ds.H500.isel(time=0).plot.pcolormesh(ax=axes[0], x="lon", y="lat", shading="flat");
ds_out_cf.isel(time=0).plot.pcolormesh(ax=axes[1], x="lon", y="lat", shading="flat");
ds_out.isel(time=0).plot.pcolormesh(ax=axes[2], x="lon", y="lat", shading="flat");
axes[0].title.set_text("Source");
axes[1].title.set_text("Target (CF)");
axes[2].title.set_text("Target");
    
for axis in axes.flatten():
    axis.coastlines();
    # zoom on north pole
    axis.set_xlim([125, 150]);
    axis.set_ylim([85, 90]);

I however encountered another problem, because the north and soutmost rows of the Merra-2 grid have only half the meridional extent, but there are no bounds provided in your input dataset. When using the conservative method to regrid, xesmf therefore uses cf_xarray to calculate/guess the bounds. The calculated bounds range from -90.25 to 90.25 degrees latitude, while the bounds of the target grid range from -90 to 90 degrees latitude. xesmf also issues a warning because of that: UserWarning: Latitude is outside of [-90, 90]

For some reason this leads to the values of the polar rows to be too high (reminding me of the problem regarding duplicated cells #109). In below example I adjusted the bounds before calculating the weights and the results were then as expected. Alternatively, one can use the conservative_normed method when calculating the weights, or the option skipna=True when applying the weights, to apply a (re)normalization.

import numpy as np
import xarray as xr
import xesmf as xe
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
ds = xr.open_dataset("grass_MERRA2_400.tavg1_2d_slv_Nx.20200101_variable.nc4")
ds = ds.cf.add_bounds("latitude")
ds = ds.cf.add_bounds("longitude")
ds["lat_bounds"] = ds["lat_bounds"].where(ds["lat_bounds"]<=90., 90.)
ds["lat_bounds"] = ds["lat_bounds"].where(ds["lat_bounds"]>=-90., -90.)

# Create weights for global target grid (in CF and non-CF representation)
regridder = xe.Regridder(ds, xe.util.grid_global(1., 1.), "conservative")
regridder_cf = xe.Regridder(ds, xe.util.grid_global(1., 1., cf=True), "conservative")

# Regrid
ds_out = regridder(ds.H500)
ds_out_cf = regridder_cf(ds.H500)

# Assert both results are equal
np.testing.assert_array_equal(ds_out.values, ds_out_cf.values)

# Plot
fig, axes = plt.subplots(nrows=3, figsize=(18,12), subplot_kw={'projection': ccrs.PlateCarree()});
ds.H500.isel(time=0).plot.pcolormesh(ax=axes[0], x="lon", y="lat", shading="flat");
ds_out_cf.isel(time=0).plot.pcolormesh(ax=axes[1], x="lon", y="lat", shading="flat");
ds_out.isel(time=0).plot.pcolormesh(ax=axes[2], x="lon", y="lat", shading="flat");
axes[0].title.set_text("Source");
axes[1].title.set_text("Target (CF)");
axes[2].title.set_text("Target");
    
for axis in axes.flatten():
    axis.coastlines();
    # zoom on north pole
    axis.set_xlim([125, 150]);
    axis.set_ylim([85, 90]);

[trimpinm]

Thank you for this response! Would the bounds still have to be specified in a similar manner when implementing bilinear/nearest neighbor regridding? Or is this method only useful when it comes to conservative regridding?

[sol1105]

For methods other than conservative*, bounds do not have to be specified, and if they are, they are ignored.