PROJ Mollweide projection singularity at the poles #1333
Description
Description
The PROJ Mollweide projection will map all points above a certain latitude (somewhere between ±89.144 and ±89.145 degrees) directly to the pole. This seems to be due to Mollweide using an iterative algorithm which is slow to converge at the poles.
As a result, this causes certain geometries to be incorrectly projected. In particular, this will affect contourf plots which use data above 89 degrees, with some polygons incorrectly being projected onto almost the entire map, blotting out previously plotted polygons.
Below is a polygon taken from a contourf. When projected, it gets 'turned inside out'.
Code to reproduce
import cartopy
import shapely.wkt
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
polygon = shapely.wkt.loads('POLYGON ((-124.6875 88.51576898014021, -124.6875 89.375, -126.5625 89.375, -128.4375 89.375, -130.3125 89.375, -132.1875 89.375, -134.0625 89.375, -135.9375 89.375, -137.8125 89.375, -139.6875 89.375, -141.5625 89.375, -143.4375 89.375, -145.3125 89.375, -147.1875 89.375, -149.0625 89.375, -150.9375 89.375, -152.8125 89.375, -154.6875 89.375, -156.5625 89.375, -158.4375 89.375, -160.3125 89.375, -162.1875 89.375, -164.0625 89.375, -165.9375 89.375, -167.8125 89.375, -169.6482699764407 89.375, -167.8125 88.48909844616483, -165.9375 88.57715321661044, -164.0625 88.66842665037734, -162.1875 88.56738315818514, -160.3125 88.75973810567709, -158.4375 88.61867874791304, -156.5625 88.48770835812856, -154.6875 88.73331726616871, -152.8125 88.6420385333048, -150.9375 89.20007552152518, -149.0625 88.46018604660821, -147.1875 88.42282488713046, -145.3125 88.39871713455433, -143.4375 88.9979823936324, -141.5625 88.60989689925717, -139.6875 88.86032315648676, -137.8125 88.58963849370639, -135.9375 88.63502159672632, -134.0625 88.53263139067644, -132.1875 88.75896594337623, -130.3125 88.76036263098899, -128.4375 88.79108697623622, -126.5625 88.52046382385734, -124.6875 88.51576898014021))')
mol = ccrs.Mollweide()
plcr = ccrs.PlateCarree()
polys = mol.project_geometry(polygon, src_crs=plcr)
plt.subplot(111)
for poly in polys:
x, y = poly.exterior.xy
plt.plot(x, y)
plt.show()