Wrong coordinates returned by `AreaDefintion.get_lonlats` for some projections in out-of-Earth locations

[ghiggi]

In the context of #546, I was testing that the new boundary extraction methodology works for all projections available in cartopy, and I stumbled in wrong values returned by AreaDefinition.get_lonlats().
After a quick investigation, the error seems to arise from wrong results by pyproj transformations !

Here below I provide a code that enable to visualize the wrong coordinates as well as snapshot of what is going on for some of the problematic projections.

Problem description

The wrong results occurs only for some CRS in out-of-Earth locations.
Usually pyproj returns np.inf values in out-of-Earth locations, while with the projections listed below, AreaDefinition.get_lonlats() returns wrong/bad longitude/latitude values. In a couple of cases, even returns latitudes outside the (-90, 90) bounds.

NOTE1: These errors causes wrong resampling with pyresample !

NOTE2: Note that for all cartopy projections not listed below, pyproj returns the correct results !

Expected Output

Usually pyproj returns np.inf values in out-of-Earth locations.

Code Sample, a minimal, complete, and verifiable piece of code

import numpy as np 
import pyproj 
import cartopy
from pyresample.future.geometry import AreaDefinition 
import matplotlib.pyplot as plt 


list_bugs = [ 
    cartopy.crs.EquidistantConic(central_longitude=0.0, central_latitude=0.0, false_easting=0.0, false_northing=0.0, standard_parallels=(20.0, 50.0), globe=None),
    cartopy.crs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=0.0, false_easting=0.0, false_northing=0.0, globe=None),
    cartopy.crs.AlbersEqualArea(central_longitude=0.0, central_latitude=0.0, false_easting=0.0, false_northing=0.0, standard_parallels=(20.0, 50.0), globe=None),
    cartopy.crs.LambertConformal(central_longitude=-96.0, central_latitude=39.0, false_easting=0.0, false_northing=0.0, standard_parallels=(33, 45), globe=None, cutoff=-30),
    cartopy.crs.EqualEarth(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.Hammer(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.Sinusoidal(central_longitude=0.0, false_easting=0.0, false_northing=0.0, globe=None),
    cartopy.crs.EckertI(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertII(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertIII(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertIV(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertV(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertVI(central_longitude=0, false_easting=None, false_northing=None, globe=None),
]

cartopy_crs = list_bugs[1]
for cartopy_crs in list_bugs:
    name = type(cartopy_crs).__name__
    x_min, x_max = cartopy_crs.x_limits
    y_min, y_max = cartopy_crs.y_limits
    area_extent = [x_min, y_min, x_max, y_max]
    proj_dict = cartopy_crs.to_dict()
    shape = (200, 300)
    area_def = AreaDefinition(
        crs=proj_dict,
        shape=shape,
        area_extent=area_extent
    )
    lons, lats = area_def.get_lonlats()
    
    plt.imshow(lons)
    plt.title(name+"Longitude"); plt.colorbar()
    plt.show()
    
    plt.imshow(lats)
    plt.title(name+" Latitude"); plt.colorbar()
    plt.show() 

To further investigate the erroneous results, just take values in the lower left corner (which are out of Earth in the specified projections) and look at the returned values:

from pyresample.geometry import get_geodetic_crs_with_no_datum_shift, TransformDirection, CRS
self = area_def
proj_wkt = self.crs_wkt

x, y = area_def.get_proj_coords(data_slice=(-1, 0))
crs = CRS.from_wkt(proj_wkt)
gcrs = get_geodetic_crs_with_no_datum_shift(crs)
transformer = pyproj.Transformer.from_crs(gcrs, crs, always_xy=True)
lon, lat = transformer.transform(x, y, direction=TransformDirection.INVERSE)
print(lon, lat)   # should be inf, inf 

Examples

Equidistant Conic

Azimuthal Equidistant

Albers Equal Area

Equal Earth

Lamber Conformal

Sinusoidal

Eckert

[djhoese]

Nice find. Some comments:

1. I think cartopy CRS objects are not pyproj subclasses so pyproj.crs.CRS.from_user_input(cartopy_crs) should work and likely uses WKT to do the conversion which is more accurate than the PROJ dict.
2. pyproj is just a wrapper around PROJ so to get this discussed/fixed upstream would require a PROJ reproducible example. Likely not that hard using the proj command line tool.
3. I don't anticipate an overwhelmingly positive response from PROJ maintainers in changing this. I semi-recently made an update so when NaNs were provided to PROJ it maintained the NaN (see LAEA projection odd handling of NaN OSGeo/PROJ#3596 and the related PR). Your case is different enough since you're passing in real values that maybe it would be something they'd consider fixing or at least merging if you did the fixes yourself.

[djhoese]

Oh another question, if you take these out of bounds lon/lats and do the inverse calculation, do you get the same expected X/Y back?

[ghiggi]

@djhoese no.
I tried: x1, y1 = transformer.transform(lon, lat, direction=TransformDirection.FORWARD) and I got different results !!!

[djhoese]

Can you paste an example here that doesn't use pyresample and shows the results of the round trip (x/y -> lon/lat -> x/y)? Preferably only pyproj in the example would be awesome.

[ghiggi]

Here it is @djhoese . A sample case for each of the problematic CRS.
I choosed testing points in the bottom left corner for all projections except for EquidistantConic, AlbersEqualArea and LambertConformal where I took an out-of-Earth disk point which is located in the middle of the first row (see above maps)

import pyproj
from pyproj.enums import TransformDirection
from pyproj import CRS

dict_problems = {'+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs': (-19970716.64831328,
  -19870474.739266966),
 '+proj=eqearth +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs': (-17186479.198676225,
  -8350962.960474116),
 '+proj=hammer +a=6378137.0 +lon_0=0 +no_defs +type=crs': (-17979962.043826804,
  -8974947.608833278),
 '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs': (-19970716.64831328,
  -9951955.900667064),
 '+proj=eck1 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs': (-18399375.367312096,
  -9184303.590539567),
 '+proj=eck2 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs': (-18399375.367312096,
  -9184303.590539567),
 '+proj=eck3 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs': (-16864798.91320002,
  -8418298.454164224),
 '+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs': (-16864798.91320002,
  -8418298.454164224),
 '+proj=eck5 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs': (-17614682.20479657,
  -8792613.107243773),
 '+proj=eck6 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs': (-17614682.20479657,
  -8792613.107243773),
 '+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=20 +lat_2=50 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs': (75949.73118667677,
  17489889.957611486),
 '+proj=aea +lat_0=0 +lon_0=0 +lat_1=20 +lat_2=50 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs': (-5841910.532119121,
  15353221.004741197),
 '+proj=lcc +lat_0=39 +lon_0=-96 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs': (58559.49249108136,
  14524620.012053927)}
 
for proj4_string, (x,y) in dict_problems.items(): 
    crs = CRS.from_proj4(proj4_string)
    gcrs = crs.geodetic_crs
    if not gcrs.prime_meridian.longitude == 0:
        raise ValueError("Just to check ...")
    transformer = pyproj.Transformer.from_crs(gcrs, crs, always_xy=True)
    lon, lat = transformer.transform(x, y, direction=TransformDirection.INVERSE)
    print(proj4_string)
    print(f"- Longitude/Latitude: ({lon},{lat})")   # should be inf, inf 
 
    x1, y1 = transformer.transform(lon, lat, direction=TransformDirection.FORWARD)
    print("- Projection coordinates:")
    print(f"  - ({x1},{y1})")
    print(f"  - ({x},{y})")

This is the result:

+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
- Longitude/Latitude: (112.9857260620742,42.71051653688908)
- Projection coordinates:
  - (8387703.848241446,8398203.2914208)
  - (-19970716.64831328,-19870474.739266966)
+proj=eqearth +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
- Longitude/Latitude: (60.075170200351245,-85.30998067351418)
- Projection coordinates:
  - (3442464.779244623,-8350962.960469991)
  - (-17186479.198676225,-8350962.960474116)
+proj=hammer +a=6378137.0 +lon_0=0 +no_defs +type=crs
- Longitude/Latitude: (14.893071680164464,-7.37221690899598)
- Projection coordinates:
  - (1646418.8895271344,-821832.8403408865)
  - (-17979962.043826804,-8974947.608833278)
+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
- Longitude/Latitude: (159.50921062391092,-89.55226021012916)
- Projection coordinates:
  - (139223.99119363955,-9951955.900667062)
  - (-19970716.64831328,-9951955.900667064)
+proj=eck1 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs
- Longitude/Latitude: (2.9850746268656287,-89.55)
- Projection coordinates:
  - (153840.93116481462,-9184303.590539567)
  - (-18399375.367312096,-9184303.590539567)
+proj=eck2 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs
- Longitude/Latitude: (2.9850746268656287,-85.31466976433765)
- Projection coordinates:
  - (153840.93116481462,-9184303.590539567)
  - (-18399375.367312096,-9184303.590539567)
+proj=eck3 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs
- Longitude/Latitude: (33.781088289342044,-89.55)
- Projection coordinates:
  - (1746407.8280480525,-8418298.454164222)
  - (-16864798.91320002,-8418298.454164224)
+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs
- Longitude/Latitude: (33.78108828934225,-85.57037018986479)
- Projection coordinates:
  - (1746407.8280480693,-8418298.45416422)
  - (-16864798.91320002,-8418298.454164224)
+proj=eck5 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs
- Longitude/Latitude: (3.9960199751023566,-89.54999999999998)
- Projection coordinates:
  - (197718.63769760213,-8792613.10724377)
  - (-17614682.20479657,-8792613.107243773)
+proj=eck6 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs
- Longitude/Latitude: (3.9960199751023566,-85.51140046588031)
- Projection coordinates:
  - (197718.63769760213,-8792613.10724377)
  - (-17614682.20479657,-8792613.107243773)
+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=20 +lat_2=50 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
- Longitude/Latitude: (-44.291137382194314,72.74673315516685)
- Projection coordinates:
  - (-1998617.4736796676,8520786.354244715)
  - (75949.73118667677,17489889.957611486)
+proj=aea +lat_0=0 +lon_0=0 +lat_1=20 +lat_2=50 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
- Longitude/Latitude: (156.8196364850849,60.133843485001925)
- Projection coordinates:
  - (6319302.333365837,12579582.00950822)
  - (-5841910.532119121,15353221.004741197)
+proj=lcc +lat_0=39 +lon_0=-96 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
- Longitude/Latitude: (-171.30536801583568,49.2011111044446)
- Projection coordinates:
  - (-4935029.573955551,3303790.1751471036)
  - (58559.49249108136,14524620.012053927)

[djhoese]

Yeah I'm having trouble understanding this. I wouldn't be surprised if the answer for some projections is "yeah, it isn't 1:1, that's just the way it is". But for the first definition in your dict I can do the inverse transform of -10-50 million in the X direction and at the reference latitude (0) and get about the same -90 longitude:

In [2]: crs = CRS.from_string('+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs')

In [4]: transformer = Transformer.from_crs(crs.geodetic_crs, crs, always_xy=True)

In [12]: transformer.transform(-10000000, 0, direction="INVERSE")
Out[12]: (-89.83152841195214, 0.0)

In [13]: transformer.transform(-50000000, 0, direction="INVERSE")
Out[13]: (-89.15764205976068, 0.0)

And somehow the -50M is closer to 0.

But then doing the lon/lat by themselves shows that it should be closer to -10 million and going all the way to -175 longitude only gives -19 million.

In [20]: transformer.transform(-89.83, 0)
Out[20]: (-9999829.857959764, 6.123129813784278e-10)

In [21]: transformer.transform(-90, 0)
Out[21]: (-10018754.171394622, 6.134717613721308e-10)

In [22]: transformer.transform(-175, 0)
Out[22]: (-19480910.888822876, 1.1928617582235876e-09)

Oh! It wraps:

In [40]: transformer.transform(-20000000, 0, direction="INVERSE")
Out[40]: (-179.66305682390427, -0.0)

In [41]: transformer.transform(-21000000, 0, direction="INVERSE")
Out[41]: (171.35379033490054, -0.0)

In [42]: transformer.transform(-22000000, 0, direction="INVERSE")
Out[42]: (162.37063749370532, -0.0)