pcolormesh with transparent colormap and Orthographic projection does not work

[AgilentGCMS]

Description

I am trying to plot a 2D field on a map with pcolormesh using a colormap with transparency. For reasons described here I am using the mplcairo backend. With that choice, plotting on (say) the Robinson projection works as expected. However, plotting on an Orthographic projection results in just the background being shown, i.e., pcolormesh has no effect.

Code to reproduce

from matplotlib import pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
import matplotlib.colors as m_colors

# create lat, lon and data arrays
lat_edges = np.linspace(-90., 90., 361)
lon_edges = np.linspace(-180., 180., 721)
lat_centers = 0.5*(lat_edges[1:]+lat_edges[:-1]) * np.pi/180.0 # convert to radian
lon_centers = 0.5*(lon_edges[1:]+lon_edges[:-1]) * np.pi/180.0 # convert to radian
lon_c, lat_c = np.meshgrid(lon_centers, lat_centers)
data = np.cos(2*lat_c) * np.sin(2*lon_c)

# create a colormap with transparency, that is blue for 0 and red for 1, with transparency in the middle
my_cmap=m_colors.LinearSegmentedColormap.from_list("", [(0,0,1,1), (1,1,1,0), (1,0,0,1)])

central_lat = 0.0
central_lon = 0.0
data_proj = ccrs.PlateCarree()

# First, orthographic projection
map_proj = ccrs.Orthographic(central_longitude=central_lon, central_latitude=central_lat)
aspect = (map_proj.x_limits[1]-map_proj.x_limits[0])/(map_proj.y_limits[1]-map_proj.y_limits[0])
height = 5.0
width  = height * aspect
fig = plt.figure(figsize=(width, height))
fig.set_facecolor('0.5')
plot_ax = plt.axes([0.,0.,1.,1.], projection=map_proj)
plot_ax.set_facecolor('k')

img = plot_ax.pcolormesh(lon_edges, lat_edges, data,
    transform=data_proj, cmap=my_cmap, vmin=-1, vmax=1, shading='flat')

# Now, Robinson
map_proj = ccrs.Robinson(central_longitude=central_lon)
aspect = (map_proj.x_limits[1]-map_proj.x_limits[0])/(map_proj.y_limits[1]-map_proj.y_limits[0])
height = 5.0
width  = height * aspect
fig = plt.figure(figsize=(width, height))
fig.set_facecolor('0.5')
plot_ax = plt.axes([0.,0.,1.,1.], projection=map_proj)
plot_ax.set_facecolor('k')

img = plot_ax.pcolormesh(lon_edges, lat_edges, data,
    transform=data_proj, cmap=my_cmap, vmin=-1, vmax=1, shading='flat')

# show all
plt.show()

Results

The Robinson projection plot shows red and blue colors at either end, with the black background showing through in the middle because the colormap is transparent there. However, the Orthographic projection just shows the black background.

Full environment definition

### Operating system Mac OSX 14.2.1 (Sonoma) ### Cartopy version 0.21.1 ### pip list ```shell Package Version -------------------------- ------------ appnope 0.1.3 asttokens 2.2.1 backcall 0.2.0 beniget 0.4.1 Brotli 1.1.0 build 1.0.3 cached-property 1.5.2 Cartopy 0.21.1 certifi 2023.11.17 cftime 1.6.3 charset-normalizer 3.3.2 contourpy 1.2.0 cppy 1.2.1 cycler 0.12.1 Cython 0.29.36 decorator 5.1.1 docutils 0.20.1 executing 1.2.0 fb-re2 1.0.7 fonttools 4.46.0 gast 0.5.2 GDAL 3.8.1 gnureadline 8.1.2 h5py 3.10.0 idna 3.4 installer 0.7.0 ipython 8.14.0 jedi 0.19.1 kiwisolver 1.4.5 lxml 4.9.3 matplotlib 3.8.2 matplotlib-inline 0.1.6 mercurial 6.6 meson 1.3.0 meson-python 0.15.0 mplcairo 0.5 mypy 1.7.1 mypy-extensions 1.0.0 netCDF4 1.6.5 numpy 1.26.1 olefile 0.46 OWSLib 0.29.3 packaging 23.2 parso 0.8.3 pexpect 4.8.0 pickleshare 0.7.5 Pillow 9.5.0 pip 23.2.1 pkgconfig 1.5.5 ply 3.11 prompt-toolkit 3.0.38 ptyprocess 0.7.0 pure-eval 0.2.2 pybind11 2.11.1 pycairo 1.23.0 pyepsg 0.4.0 Pygments 2.15.1 pyparsing 3.1.1 pyproj 3.6.0 pyproject_hooks 1.0.0 pyproject-metadata 0.7.1 pyshp 2.1.3 python-dateutil 2.8.2 pythran 0.14.0 pytz 2023.3.post1 PyYAML 6.0.1 requests 2.31.0 roman 4.1 SciPy 1.10.1 setuptools 68.2.2 setuptools-scm 8.0.4 setuptools-scm-git-archive 1.1 shapely 2.0.1 six 1.16.0 stack-data 0.6.2 tqdm 4.66.0 traitlets 5.9.0 types-psutil 5.9.5.16 types-setuptools 68.2.0.2 typing_extensions 4.9.0 unicodedata2 15.1.0 urllib3 1.26.18 wcwidth 0.2.12 wheel 0.42.0 zopfli 0.2.3 ```

[AgilentGCMS]

Did you try with the mplcairo backend? For reasons mentioned in the ticket, I have to use the mplcairo backend, because otherwise I get a moire pattern.

On the mplcairo backend, Orthographic and Robinson have different behaviors,

On Thu, Dec 21, 2023 at 11:11 AM Ruth Comer @.***> wrote:

I just tried the above code with Cartopy 0.22 on Linux, and I got this for the Orthographic projection:

image.png (view on web) https://github.com/SciTools/cartopy/assets/10599679/06c9068d-6a85-4322-995a-745446d99629

— Reply to this email directly, view it on GitHub https://github.com/SciTools/cartopy/issues/2305#issuecomment-1866581485, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA5QCQL3ULCS5CXXYXTV3NLYKRNUZAVCNFSM6AAAAABA6S5AEKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRWGU4DCNBYGU . You are receiving this because you authored the thread.Message ID: @.***>

[rcomer]

Sorry, I deleted my comment because I realised too late that the backend was probably the difference.

[greglucas]

I wonder if this has to do with the transform of points being nan's... map_proj.transform_point(120, 20, data_proj) == (nan, nan) @anntzer do you know how mplcairo handles nan output from a non-affine transform? Perhaps Agg handles the nan's in the quadmesh paths differently than mplcairo.

Changing: lon_edges = np.linspace(-180., 180., 721) to the valid transform range for longitudes: np.linspace(-90, 90, 721) gives a good plot.

[anntzer]

I didn't really check how things get threaded to the backend, but mplcairo doesn't do any specific checking for nans in draw_quad_mesh and it seems likely that passing a nan will just break cairo. It would be nice if the backend spec (essentially the docstrings of the base methods in RendererBase) specified clearly which methods need to support nans, and for which inputs. I assume here it's for the coordinates parameter? Adding a nan check seems doable, something like (untested)

diff --git i/ext/_mplcairo.cpp w/ext/_mplcairo.cpp
index 2b46997..dbca760 100644
--- i/ext/_mplcairo.cpp
+++ w/ext/_mplcairo.cpp
@@ -1506,6 +1506,16 @@ void GraphicsContextRenderer::draw_quad_mesh(
     auto const& pattern = cairo_pattern_create_mesh();
     for (auto i = 0; i < mesh_height; ++i) {
       for (auto j = 0; j < mesh_width; ++j) {
+        if (!(std::isfinite(coords_raw(i, j, 0))
+              && std::isfinite(coords_raw(i, j, 1))
+              && std::isfinite(coords_raw(i, j + 1, 0))
+              && std::isfinite(coords_raw(i, j + 1, 1))
+              && std::isfinite(coords_raw(i + 1, j, 0))
+              && std::isfinite(coords_raw(i + 1, j, 1))
+              && std::isfinite(coords_raw(i + 1, j + 1, 0))
+              && std::isfinite(coords_raw(i + 1, j + 1, 1)))) {
+          continue;
+        }
         cairo_mesh_pattern_begin_patch(pattern);
         cairo_mesh_pattern_move_to(
           pattern, coords_raw(i, j, 0), coords_raw(i, j, 1));

(and the same in the draw-one-at-a-time case just above).

[AgilentGCMS]

I can confirm that this only occurs if there are lat/lon points that are not shown on the plot due to the projection, as @greglucas suggested. Even a complicated projection like InterruptedGoodeHomolosine displays fine, while any projection that hides part of the globe, like NearsidePerspective, creates pure black images. In the above example, changing the longitudes to only those that can be displayed, by lon_edges = np.linspace(-90., 90., 361), fixes the problem.

One additional diagnostic is that even one or two longitudes outside the displayed area renders the entire image black. So e.g., lon_edges = np.linspace(-90.5, 90., 362), which puts just one longitude outside the display area, creates a pure black image with the Orthographic projection.

[greglucas]

This is fixed and merged in upstream in mplcairo now. Closing here as there is nothing to do within Cartopy.