Skip to content

Possible regression: contourf() with 1D X/Y flips vertically in Cartopy 0.25.0 when masked data present at edge #2620

@tomvothecoder

Description

@tomvothecoder

Overview

Hi Cartopy team, thanks for maintaining Cartopy. While testing downstream diagnostics in e3sm_diags, we noticed a regression in contourf() behavior introduced in Cartopy 0.25.0. This issue is adapted from the original e3sm_diags report (here) and includes a minimal reproducible example and a narrowed commit window to help investigate.

Cartopy version: 0.25.0
Last known good: 0.24.0
Introduced by: eb4d042
Matplotlib: 3.10.8
Platform: Linux / macOS (reproducible)

What happened?

For our test case, Cartopy 0.25.0 and contourf() produces a vertically flipped (Y-axis–inverted) result when:

  • X and Y are provided as 1D coordinate arrays
  • Data is 2D (y, x)
  • Masked values are present along one edge (e.g., southernmost latitudes)
  • A geographic transform is used (e.g., PlateCarree)

This behavior does not occur in Cartopy 0.24.0 and earlier. The issue appears most clearly in difference plots (e.g., test - ref), where masked regions are common.

I think the output of contourf() using 1D X/Y coordinates should remain consistent between Cartopy 0.24.x and 0.25.x.

Regression window

The issue first appears after this commit:

The prior commit works correctly:


Minimal Reproducible Example (MVCE)

Note: The data includes masked values at the grid edge.
This script below was adapted based on how e3sm_diags produces plots. Apologies if it is a bit verbose.

Details
import urllib.request
import os

import cartopy
import matplotlib
import numpy as np
import xarray as xr
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter

matplotlib.use("Agg")
import matplotlib.pyplot as plt  # isort:skip  # noqa: E402

# Download the data
# --------------------------------------------------------------------------
# Source: /lcrc/group/e3sm/public_html/diagnostic_output/ac.tvo/tests/1026_cartopy_0250_mvce/
public_path = "https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.tvo/tests/1026_cartopy_0250_mvce/"
local_data_dir = "./data_cartopy_0250_mvce"
os.makedirs(local_data_dir, exist_ok=True)

def download_data():
    files = ["x.nc", "y.nc", "var.nc"]

    for fname in files:
        url = public_path + fname
        local_path = os.path.join(local_data_dir, fname)

        if not os.path.exists(local_path):
            print(f"Downloading {url} to {local_path} ...")
            # Set a User-Agent header to avoid HTTP 403 errors
            req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0'})
            with urllib.request.urlopen(req) as response, open(local_path, 'wb') as out_file:
                out_file.write(response.read())
        else:
            print(f"File {local_path} already exists.")

download_data()

# Open the data
# --------------------------------------------------------------------------
x = xr.open_dataarray(os.path.join(local_data_dir, "x.nc"))
y = xr.open_dataarray(os.path.join(local_data_dir, "y.nc"))
var = xr.open_dataarray(os.path.join(local_data_dir, "var.nc"))

# Create the figure and axis with Cartopy projection
# --------------------------------------------------------------------------
fig = plt.figure(figsize=(8.5, 11.0), dpi=150)
projection = cartopy.crs.PlateCarree(central_longitude=180)
ax = fig.add_axes((0.1691, 0.1112, 0.6465, 0.2258), projection=projection)

# Set the longitude and latitude limits.
lon_west = 0
lon_east = 360
lat_south = -90
lat_north = 90

ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection)

# Create the contour plot
# --------------------------------------------------------------------------
transform = cartopy.crs.PlateCarree()
contour_plot = ax.contourf(
    x,
    y,
    var,
    cmap="BrBG",
    transform=transform,
    norm=None,
    levels=None,
    extend="both",
)

# Configure aspect ratio and coast lines.
# --------------------------------------------------------------------------
ax.set_aspect((lon_east - lon_west) / (2 * (lat_north - lat_south)))
ax.coastlines(lw=0.3)

# Configure x and y axes:
# --------------------------------------------------------------------------
x_ticks = np.array([0., 60., 120., 180., 240., 300., 359.5])
y_ticks = np.array([-90, -60, -30, 0, 30, 60, 90])
ax.set_xticks(x_ticks, crs=transform)
ax.set_yticks(y_ticks, crs=transform)
ax.tick_params(labelsize=8.0, direction="out", width=1)
ax.xaxis.set_ticks_position("bottom")
ax.yaxis.set_ticks_position("left")

# Add longitude and latitude formatters
lon_formatter = LongitudeFormatter(
    zero_direction_label=True, number_format=".0f"
)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

# Add color bar
# --------------------------------------------------------------------------
cbax_rect = (0.8326, 0.13269999999999998, 0.0326, 0.1792)
cbax = fig.add_axes(cbax_rect)
cbar = fig.colorbar(contour_plot, cax=cbax)
cbar.ax.tick_params(labelsize=9.0, length=0)

# Save the figure
# --------------------------------------------------------------------------
conda_env = os.environ.get("CONDA_DEFAULT_ENV", "unknown_env")
output_path = os.path.abspath(f"{conda_env}_mvce.png")
fig.savefig(output_path)

print(f"Saved figure at {output_path}")

Observed workaround (not ideal)

Using a 2D meshgrid and transform_first=True avoids the inversion:

x2d, y2d = np.meshgrid(x, y)
ax.contourf(
    x2d,
    y2d,
    var,
    transform=cartopy.crs.PlateCarree(),
    transform_first=True,
)

However, this workaround produces unexpected results for other plots in our diagnostic sets.

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions