Mollweide projection - completely corrupted contourf when there are gridpoints at the pole

[jacopok]

Description

This is a similar issue to #1333 or #2016, but the effect is much more dramatic, not only limited to the poles.
It seems to happen when there are gridpoints exactly at the poles, not when they are further north than 89.14,
so #1334 might not fix it.
Therefore, I'm thinking that it might be a different problem.

Basically, plotting the exact same data, this is what I'm getting:

Setting the latitude grid to go from -89.9 to 89.9 degrees fixes the problem, but it still seems like a bug.

Code to reproduce

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

def function(theta, phi):
    return 0.5 * (1 + np.cos(theta) ** 2) * np.cos(2 * phi)

nlon = 240
nlat = 120
lon_degrees = np.linspace(0, 360.0, nlon)
lat_degrees = np.linspace(-90., 90., nlat)

lon = np.deg2rad(lon_degrees)
lat = np.deg2rad(lat_degrees)

lon2d, lat2d = np.meshgrid(lon, lat)

data = function(np.pi/2 - lat2d, lon2d)

fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(2, 2, 1)

mappable = ax1.contourf(
    np.rad2deg(lon2d),
    np.rad2deg(lat2d),
    data,
    cmap=plt.get_cmap("seismic"),
    levels=np.linspace(-1, 1, num=50),
)
plt.colorbar(mappable, ax=ax1)
ax1.set_title('Bare contourf')

proj = ccrs.Mollweide()
ax2 = fig.add_subplot(2, 2, 2, projection=proj)
mappable = ax2.contourf(
    np.rad2deg(lon2d),
    np.rad2deg(lat2d),
    data,
    cmap=plt.get_cmap("seismic"),
    levels=np.linspace(-1, 1, num=50),
    transform=ccrs.PlateCarree(),
)
ax2.coastlines()
ax2.set_global()
plt.colorbar(mappable, ax=ax2)
ax2.set_title('Mollweide')

proj = ccrs.PlateCarree()
ax3 = fig.add_subplot(2, 2, 3, projection=proj)
mappable = ax3.contourf(
    np.rad2deg(lon2d),
    np.rad2deg(lat2d),
    data,
    cmap=plt.get_cmap("seismic"),
    levels=np.linspace(-1, 1, num=50),
    transform=ccrs.PlateCarree(),
)
ax3.coastlines()
ax3.set_global()
plt.colorbar(mappable, ax=ax3)
ax3.set_title('PlateCarree')

proj = ccrs.Robinson()
ax4 = fig.add_subplot(2, 2, 4, projection=proj)
mappable = ax4.contourf(
    np.rad2deg(lon2d),
    np.rad2deg(lat2d),
    data,
    cmap=plt.get_cmap("seismic"),
    levels=np.linspace(-1, 1, num=50),
    transform=ccrs.PlateCarree(),
)
ax4.coastlines()
ax4.set_global()
plt.colorbar(mappable, ax=ax4)
ax4.set_title('Robinson')

Versions

Cartopy version 0.21.0, matplotlib version 3.6.1.

[lgolston]

Thanks for the report - a bit hard to figure out what is going on here. I have not found the issue but did notice that the Mollweide plot generally works here without adjusting the latitude if using a smaller numbers of color levels (up to num = 25, but also num=30). I wonder if that points to what is going on? It also works for any number of levels with transform_first=True.

[lgolston]

After seeing #1076, this problem appears to be similar. Was able to create some test cases based on your example to better show what is going on. Each case isolates a specific path which causes issues (there could be multiple). It also shows a comparison to the same plot, but with the endpoints of the lat grid shifted slightly away from the pole.

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib

proj = ccrs.Mollweide()

def function(theta, phi):
    return 0.5 * (1 + np.cos(theta) ** 2) * np.cos(2 * phi)

def plot_Mollweide(ax,dlat,num_levels):
    nlon = 240
    nlat = 120
    lon_degrees = np.linspace(0, 360.0, nlon)
    lat_degrees = np.linspace(-90, 90, nlat)
    lat_degrees[0] += dlat
    lat_degrees[-1] -= dlat
    
    lon = np.deg2rad(lon_degrees)
    lat = np.deg2rad(lat_degrees)
    
    lon2d, lat2d = np.meshgrid(lon, lat)
    
    data = function(np.pi/2 - lat2d, lon2d)
    
    levels = np.linspace(-1,1,num=num_levels)

    mappable = ax.contourf(
        np.rad2deg(lon2d),
        np.rad2deg(lat2d),
        data,
        cmap=plt.get_cmap("seismic"),
        levels=levels,
        transform=ccrs.PlateCarree()
    )
    return mappable

def get_path(mappable,path_ii,path_jj):
    # select specific path
    for ii, pc in enumerate(mappable.collections):
        if ii==path_ii:
            for jj, path in enumerate(pc.get_paths()):     
                if jj == path_jj:
                    return path

# %% create test cases
test_case = 3

if test_case == 1:
    dlat = 1E-4
    num_levels = 18
    path_ii = 10
    path_jj = 1
elif test_case == 2:
    dlat = 1E-4
    num_levels = 28
    path_ii = 20
    path_jj = 1
elif test_case == 3:
    dlat = 1E-6
    num_levels = 30
    path_ii = 23
    path_jj = 0 # can be paths 0 and 5 for many values of ii here

# %% compare two examples, only adjusting the latitude endpoint slightly (dlat)
# create figure
fig = plt.figure(figsize=(8, 5))
plot_proj = ccrs.Mollweide(central_longitude=0)

# contourf plots
ax1 = fig.add_subplot(2, 2, 1, projection=plot_proj)
c1 = plot_Mollweide(ax1, 0, num_levels)

ax2 = fig.add_subplot(2, 2, 2, projection=plot_proj)
c2 = plot_Mollweide(ax2, dlat, num_levels)

# extract a specific path from the collection
path_c1 = get_path(c1, path_ii, path_jj)
path_c2 = get_path(c2, path_ii, path_jj)

# plot as PathPatch
ax3 = fig.add_subplot(2, 2, 3, projection=plot_proj)
ax4 = fig.add_subplot(2, 2, 4, projection=plot_proj)
ax3.add_patch(matplotlib.patches.PathPatch(path_c1, transform=ccrs.PlateCarree()))
ax3.coastlines()
ax3.set_global()
ax4.add_patch(matplotlib.patches.PathPatch(path_c2, transform=ccrs.PlateCarree()))
ax4.coastlines()
ax4.set_global()
fig.suptitle('test case: ' + str(test_case) + '; num_levels = ' + str(num_levels))

[lgolston]

Working through the plotting functions, the problem arises from _project_linear_ring in crs.py and ultimately project_linear in trace.pyx. Here is a more reduced example, which replicates the behavior of test case 1:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.trace
import shapely.geometry as sgeom
import matplotlib

# create test data
max_lat = 90
lon1 = 140
lon2 = 150
N = 120
x1 = np.linspace(lon1, lon1, num=N)
x2 = np.linspace(lon2, lon2, num=N)
y1 = np.linspace(-max_lat,max_lat,num=N)
y2 = np.linspace(max_lat,-max_lat,num=N)

vertices = np.vstack( ([lon1, -max_lat],
                       np.vstack((x1,y1)).T,
                       np.vstack((x2,y2)).T,
                       [lon1, -max_lat]) )

# transform and plot
src_crs = ccrs.PlateCarree()
dest_crs = ccrs.Mollweide()
plot_proj = ccrs.Mollweide(central_longitude=0)

path = matplotlib.path.Path(vertices,closed=True)

fig = plt.figure(figsize=(6, 2))

# left subplot: with patch
ax1 = fig.add_subplot(1, 2, 1, projection=plot_proj)
ax1.add_patch(matplotlib.patches.PathPatch(path, transform=src_crs))
ax1.set_global()
ax1.set_title('Patch')

# right subplot: with low-level plotting
ax2 = fig.add_subplot(1, 2, 2, projection=ccrs.Mollweide(central_longitude=0))

linear_ring = sgeom.LinearRing(vertices)
multi_line_string = cartopy.trace.project_linear(linear_ring, src_crs, dest_crs)

colors = ['blue','red','green','magenta','black','yellow']
for count, ls in enumerate(multi_line_string.geoms):
    ax2.add_geometries(ls, dest_crs, color=colors[count], alpha=1)
    print(ls.is_closed)

ax2.set_global()
ax2.set_title('multi_line_string')

fig.suptitle('N=120')
fig.savefig('reduced_example.png')

With N = 100 or less, project_linear returns a single closed geometry.
With N = 120, there are two non-closed geometries, but _project_linear_ring handles it to keep the patch from filling the globe
With N = 140 and beyond, there are four non-closed geometries and the patch fills the globe

There is an open pull request #1739 that I wonder if would help, but there @kdpenner mentioned suspecting an issue is in trace which is also what I am seeing here.

[dopplershift]

This seems related to the behavior in #2240 around the bounds/boundary for Geostationary. If I increased the number of points, I started getting errors. I agree something isn't quit right in trace.project_linear(), but I'm without cycles to try to dig in.

[kdpenner]

Thanks for reminding me about my PR. I didn't expect it to go anywhere because I felt like I wasn't addressing the root issues but I am more than happy to implement suggested changes.

[lgolston]

Okay, I was able to track this down through trace (notably seeing POINT_OUT states in bisect once the problem starts happening) and ultimately back up to the definition of the boundary in _WarpedRectangularProjection in crs. I did not realize the boundary is discrete, so you can have points along the paths generated by countorf that register as being outside of the domain, e.g.:

import cartopy.crs as ccrs
from shapely.geometry import Point
dest_crs = ccrs.Mollweide()
test_point = dest_crs.transform_point(100, -89.4, ccrs.PlateCarree())
print(dest_crs.domain.contains(Point(test_point))) # returns False

Simply bumping n = 91 to n = 121 at

cartopy/lib/cartopy/crs.py

Lines 2095 to 2106 in d1b43fb

	# Obtain boundary points
	minlon, maxlon = self._determine_longitude_bounds(central_longitude)
	n = 91
	lon = np.empty(2 * n + 1)
	lat = np.empty(2 * n + 1)
	lon[:n] = minlon
	lat[:n] = np.linspace(-90, 90, n)
	lon[n:2 * n] = maxlon
	lat[n:2 * n] = np.linspace(90, -90, n)
	lon[-1] = minlon
	lat[-1] = -90
	points = self.transform_points(self.as_geodetic(), lon, lat)

fixes the original plot from @jacopok. Still need to investigate how best to implement this (e.g. maybe just add extra points near the poles for Mollweide? Do other projections need this?).