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. https://github.com/SciTools/cartopy/issues/2336#issuecomment-1978870581), 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 ;)