Wrong contour filling when rotating Stereographic projection

[guidocioni]

I'm really sorry if this is not a MWE but it is not easy to try and recreate a similar dataset, and I've spent already few hours trying to narrow down this issue :) The code used to make the plot is at the end of this post.

I'm plotting data that is natively defined in a Stereographic projection with a central_longitude=-80. In order to make a plot over Europe I'm rotating it to central_longitude=10.

This is the resulting plot when using the original projection with central_longitude=-80

This is the resulting plot when using central_longitude = 10 but leaving everything else the same

Notice the completely different (wrong) colors.

I'm not sure if there's something wrong going on with the transformation of polygons, or whether this has to do with the masked array, or the stereographic projection.

Here is the minimal code needed to reproduce this problem. Again I'm sorry it's not self-contained but I can upload the data if needed.

fig, ax = plt.subplots(
    figsize=(12, 12),
    subplot_kw={
        "projection": ccrs.Stereographic(
            central_latitude=90, central_longitude=-80, true_scale_latitude=60
        )
    },
)
extent_lonlat = [-17, 41, 35, 72]
scale='50m'
ax.set_extent(extent_lonlat, crs=ccrs.PlateCarree())
ax.add_feature(cfeature.OCEAN.with_scale(scale), facecolor="#82aac7")
ax.add_feature(
    cfeature.LAND.with_scale(scale),
    facecolor=(0.8, 0.7137254901960784, 0.615686274509804),
)

x = subset_plot["x"].values
y = subset_plot["y"].values
X, Y = np.meshgrid(x, y)
data = subset_plot.where(subset_plot > 5).values

cf = ax.contourf(
    X,
    Y,
    data,
    cmap=cmap,
    norm=norm,
    levels=levels_snow,
    transform=dataset_crs,
)

[guidocioni]

Here is an animation of a zoomed-in region where the central_longitude angle is rotated from -80 degrees to 70 degrees. Somewhere between -10 and 0 degrees one of the color of the contours fill most of the plot covering land, ocean and parts of the contours.

[rcomer]

Hi @guidocioni, what Cartopy version are you using? We did fix some of these issues at v0.25 but unfortunately there are plenty remaining.

[guidocioni]

Hi @rcomer , I'd have to check because I don't remember anymore on which device I was working when producing this particular map, but I'm pretty certain it is a fairly recent version; as far as I can see 0.25 was released back in August, right? So for sure I'm using that version.

Is there any workaround available to force the right behaviour? Right now I'm stuck with pcolormesh but that's less than optimal for some use cases...

[rcomer]

Hi @guidocioni the problem can be very sensitive to the exact numbers involved (see e.g. #2336 (comment)), so you may find that adding a tiny number to your levels_snow makes a difference.

[guidocioni]

Hi @guidocioni the problem can be very sensitive to the exact numbers involved (see e.g. #2336 (comment)), so you may find that adding a tiny number to your levels_snow makes a difference.

ok well, I'll try a montecarlo simulation to find the optimal level values to not cause the filling issue 😁

[greglucas]

Have you tried the transform_first=True keyword argument?

https://cartopy.readthedocs.io/stable/gallery/scalar_data/contour_transforms.html#sphx-glr-gallery-scalar-data-contour-transforms-py

[guidocioni]

Have you tried the transform_first=True keyword argument?

https://cartopy.readthedocs.io/stable/gallery/scalar_data/contour_transforms.html#sphx-glr-gallery-scalar-data-contour-transforms-py

Hey @greglucas , good idea, I didn't think about trying out that one. I can say it is indeed much better, although there are still some angles (e.g. central_longitude=10) that trigger weird behaviours

But at least now I can get a nice plot for most of the other angles that didn't work before!
Here the same image for central_longitude=9

BTW, as I'm usually plotting directly from xarray by using the .plot accessor, do you know if there's a way to pass that keyword to the contourf call? Right now if I try to do

da.plot.contourf(
    ax=ax,
    transform=dataset_crs,
    transform_first=True,
)

I get

ValueError: The X and Y arguments must be gridded 2-dimensional arrays

but if I try without transform_first or just using ax.contourf then everything works fine.
I understand this may be out of the scope of cartopy but maybe you already know a workaround ;)

[tomvothecoder]

Hey there — I’m wondering if this might be related to an issue I’m seeing in my repository (thread).

I compute a difference field between test and reference data, and the resulting array includes np.nan values along the bottom edge (masked).

• With Cartopy v0.24.0, the third plot renders with correct colors (plot link).
• With Cartopy v0.25.0, the same plot appears inverted in latitude (plot link).
  • If I flip the latitude coordinates (descending) in Xarray, the plot displays correctly.

Investigation

After some debugging, I narrowed the issue down to calling contourf() with 1-D X/Y coordinates + masked data. It seems this may stem from behavior introduced in SciTools/cartopy#2461.

If I instead pass 2-D X/Y coordinates to contourf() and set transform_first=True, the inversion disappears (plot link). However, we're concerned about the side effects of transform_first=True, since the docs warn:

"This does have a negative impact on the wrapped coordinates as one can see in the second axes that the data does not extend to the full global extent."

So I’m wondering: Is this a bug, and is it related to the issue discussed here?

My current hypothesis:

When using contourf() with 1-D x/y coordinates and masked values (such as
diff fields), Cartopy’s default transform path can misinterpret the grid
orientation and produce an inverted-looking plot—even though the latitude
coordinates and data are correct.

Converting x/y to explicit 2-D meshgrids and enabling `transform_first=True`
forces contourf() to use Cartopy’s polygon-based transformation path instead
of the raster/image path affected by the half-pixel shift. This restores the
correct visual orientation for diff plots while retaining the improved
top/bottom handling introduced in Cartopy 0.25.

[greglucas]

@tomvothecoder, I'm not sure I'm fully following your post because there are a few factors at play there.

I don't understand how the half-pixel shift would impact this because that is only used by the image sampling code, so I don't think that is the culprit (were you somehow able to bisect to that commit?). All contour paths are projected using the same code path, the difference with transform_first=True is in when the contouring is done. If you are going from a lat/lon data grid to say a Robinson map if transform_first=False the contouring is done in lat/lon (data) space and then those contour lines are projected into the Robinson map, if transform_first=True the individual data points are first projected into the Robinson map and then the contouring algorithm takes place producing lines directly in the Robinson projection.

[tomvothecoder]

Hi @greglucas, I appreciate your prompt and informative response. I will do more debugging and try to produce a minimum reproducible example.

I don't understand how the half-pixel shift would impact this because that is only used by the image sampling code, so I don't think that is the culprit (were you somehow able to bisect to that commit?).

The half-pixel shift PR actually seems to have fixed some of our other test and reference data plots with Y-coordinate masking at the bottom and top of images. I guessed that it might also have affected the difference plot by inverting plot colors on the Y-axis some how.

It is a great idea to actually verify my guess by bisecting at that specific commit. Will do.

All contour paths are projected using the same code path, the difference with transform_first=True is in when the contouring is done. If you are going from a lat/lon data grid to say a Robinson map if transform_first=False the contouring is done in lat/lon (data) space and then those contour lines are projected into the Robinson map, if transform_first=True the individual data points are first projected into the Robinson map and then the contouring algorithm takes place producing lines directly in the Robinson projection.

Ah, thanks for explaining what transform_first actually does. This makes sense to me now.