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

[ghiggi]

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

[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 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)