Inconsistent transformations between data and features with PROJ 9.8

[aulemahal]

Probably related to #2634.

I went on Google Maps and took the longitude and latitude of a spot near the point of an island near where I live. The exact spot isn't important , I just wanted to have something we can easily find on a map.

Then I plotted an area around that spot, adding the features for land and coastlines. The first plot is done in a projection of LambertConformal using a spherical globe and the second one is done with LambertConformal but with the default globe, which has a different radius than the previous one.

My issue is that the result is different.

import cartopy.crs as ccrs
import cartopy.feature as cfeat
import matplotlib.pyplot as plt

sph = ccrs.Globe(ellipse='sphere', semimajor_axis=6370997)
# The exact same projection as the Land and Coastline data, this is taken from their `*.prj` file.
LL = ccrs.Projection('GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]')
# Spherical versions of the projections
LCs = ccrs.LambertConformal(globe=glb)
PCs = ccrs.PlateCarree(globe=glb)
# Ellipsoidal versions
PCe = ccrs.PlateCarree()
LCe = ccrs.LambertConformal(globe=PCe.globe)

Map 1 : Spherical version of Lambert Conformal

fig, ax = plt.subplots(subplot_kw={'projection': LCs})
ax.set_extent((1536000, 1973000, 856000, 1249000), crs=LCs)

# The actual point
ax.plot(-71.1021009, 46.8643316, 'ro', transform=LL)
# What if I cheat ?
ax.plot(-71.1021009, 46.8643316, 'bo', transform=PCe)

ax.add_feature(cfeat.LAND, color='gray')
ax.coastlines()

Map 2 : Ellipsoidal version of Lambert Conformal

fig, ax = plt.subplots(subplot_kw={'projection': LCe})
ax.set_extent((1536000, 1973000, 856000, 1249000), crs=LCe)

# The actual point
ax.plot(-71.1021009, 46.8643316, 'ro', transform=LL)
# What if I cheat ?
ax.plot(-71.1021009, 46.8643316, 'bo', transform=PCe)

ax.add_feature(cfeat.LAND, color='gray')
ax.coastlines()

I think map 1 is incorrect. Assuming Google Maps gives me a geographical coordinate and not a projected one, I would assume the transformation to work in the same way for my point than for the features. Here it appears to be done differently ?

The blue point is if I assume the input coordinates are actually in the EQC projection with an ellipsoid. This is the new behaviour that comes with PROJ 9.8. Previously, the default PlateCarree was an EQC projection but with a spherical earth. The blue point should be wrong in both maps, no ?