Open ghiggi opened 10 months ago
Nice find. Some comments:
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.proj
command line tool.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?
@djhoese no.
I tried: x1, y1 = transformer.transform(lon, lat, direction=TransformDirection.FORWARD)
and I got different results !!!
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.
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)
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)
In the context of https://github.com/pytroll/pyresample/pull/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 bypyproj
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
returnsnp.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
returnsnp.inf
values in out-of-Earth locations.Code Sample, a minimal, complete, and verifiable piece of code
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:
Examples
Equidistant Conic
Azimuthal Equidistant
Albers Equal Area
Equal Earth
Lamber Conformal
Sinusoidal
Eckert