Inconsistent plots between uniform grid and non-uniform climate model grid

[ayghri]

Description

I'm using the default example with PlateCarree to plot the SSH map form CESM2 model using non-uniformly spaced longitudes/latitudes. The plots are different from what I get when I run a interpolation to 1 °-grid then plot the map. I used the same levels.

Is this a feature or a bug?

The non-uniform grid:

The contour filling using non-uniform grid:

The plot after interpolating to uniform 1°-grid:

Code to reproduce

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy import feature, util
import numpy as np
from netCDF4 import Dataset
from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator

# SSH from CESM2
ssh_data = Dataset(
    "/simulations/cesm2/SSH/b.e21.BHISTsmbb.f09_g17.LE2-1011.001.pop.h.SSH.195001-195912.nc"
)
# POP model to extract unmasked lat/long
pop1d = Dataset("/data/pop/map_1x1d_to_gx1v6_bilin_da_100716.nc")
tlongs = pop1d.variables["xc_b"][:].reshape(384, 320)
tlats = pop1d.variables["yc_b"][:].reshape(384, 320)
values = ssh_data.variables["SSH"][:][0].data
mask = ssh_data.variables["SSH"][:][0].mask
tgrid = np.stack([tlongs, tlats]).transpose(1, 2, 0)[np.logical_not(mask)]

# parameters of the plot
resolution = "110"
levels = np.linspace(-250, 200, 26)

# plot using example with original CESM2 grid
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=210))
plt.title("SSH Plot of Jan 1950 from CESM2 using original non-uniform grid")
plt.contourf(
    tlongs,
    tlats,
    np.ma.masked_where(mask, values),
    levels=levels,
    transform=ccrs.PlateCarree(),
    corner_mask=True,
    cmap="coolwarm",
    antialiased=False,
)
ax.add_feature(feature.LAND)
ax.coastlines(resolution=f"{resolution}m", linewidth=1)
ax.gridlines(draw_labels=True, linestyle="--", color="black")
cbar = plt.colorbar(shrink=0.5)
plt.show()

# Changes CESM2 [0,360] longs to [-180,180]
tlongs_180 = tlongs - 360 + (tlongs < 180) * 360
# Create uniform grid
uniform_long = np.arange(-360, 360, 2) / 2
uniform_lat = np.arange(-180, 180, 2) / 2
uniform_longs, uniform_lats = np.meshgrid(uniform_long, uniform_lat)
# Linear interpolation part
padded_tlongs = np.concatenate([tlongs_180 - 360, tlongs_180, tlongs_180 + 360], axis=1)
padded_tlats = np.concatenate([tlats, tlats, tlats], axis=1)
padded_map = np.concatenate([values, values, values], axis=1)
padded_mask = np.concatenate([mask, mask, mask], axis=1) * 1.0
coords = np.stack([padded_tlongs, padded_tlats]).transpose(1, 2, 0).reshape(-1, 2)
uni_coords = np.stack([uniform_longs, uniform_lats]).transpose(1, 2, 0).reshape(-1, 2)
## interpolating the data
map_interpolator = LinearNDInterpolator(coords, padded_map.reshape(-1))
uniform_map = map_interpolator(uni_coords)
uniform_map = uniform_map.reshape(uniform_longs.shape[0], uniform_lats.shape[1])
## interpolating the mask
mask_interpolator = LinearNDInterpolator(coords, padded_mask.reshape(-1))
uniform_mask = mask_interpolator(uni_coords)
uniform_mask = uniform_mask.reshape(uniform_longs.shape[0], uniform_lats.shape[1])
uniform_mask = np.logical_not(uniform_mask == 0)

# plotting the interpolated map
uniform_map = np.ma.masked_where(uniform_mask, uniform_map)
cyclic_map, cyclic_long = util.add_cyclic_point(uniform_map, uniform_long)
## plot using example
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=210))
plt.title("SSH Plot of Jan 1950 from CESM2 using 1d-grid interpolation")
plt.contourf(
    cyclic_long,
    uniform_lat,
    cyclic_map,
    levels=levels,
    transform=ccrs.PlateCarree(),
    corner_mask=True,
    cmap="coolwarm",
    antialiased=False,
)
ax.add_feature(feature.LAND, facecolor="white")
ax.coastlines(resolution=f"{resolution}m", linewidth=1)
ax.gridlines(draw_labels=True, linestyle="--", color="black")
cbar = plt.colorbar(shrink=0.5)
plt.show()

Full environment definition

Operating system

Linux

Cartopy version

pip list

python==3.11
Cartopy==0.22.0
contourpy==1.2.0

[dopplershift]

Looks like a bug to me. Tracking this down could really use a simplified case or direct access to the data to understand what's going on.

[rcomer]

I'm going to close this one as we cannot run the example without the data and therefore cannot investigate.