Comments (6)
Nice find. Some comments:
- 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. - 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. - 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.
from pyresample.
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?
from pyresample.
@djhoese no.
I tried: x1, y1 = transformer.transform(lon, lat, direction=TransformDirection.FORWARD)
and I got different results !!!
from pyresample.
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.
from pyresample.
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)
from pyresample.
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)
from pyresample.
Related Issues (20)
- `get_neighbour_info` slows down significantly when working with large target rasters using many segments HOT 3
- `gradient_search` fails when resampling Himawari data HOT 7
- EWA resampling in 1.27 slows down four times than 1.26.1 HOT 69
- Refactor boundary creation logic
- area definition for a rotated pole coordinate system HOT 24
- Index Error when calling `boundary` with non full disk geos ara HOT 5
- Bug where boundary coordinates wrap around 180 degrees when geostationary AreaDefinition extends beyond 180 degrees HOT 7
- Upgrade to Cython 3.0 and check annotations
- Compatibility with libproj v9.3 HOT 23
- How should this warning be addressed? HOT 2
- Catch or fix pyproj UserWarning when loading an AreaDefinition from a netCDF/CF file HOT 1
- Make a Swath definition subclass for interpolated data HOT 4
- The finding of the optimal radius of influence makes assumption on the ordering of the dimensions in the longitude array HOT 4
- dump DynamicAreaDefinition to YAML
- errors in area definition should not be silently ignored HOT 1
- Resampling GOES mesoscale data to my area gives blank data HOT 10
- Remove Configobj as Dependency Due To Security Vulnerability and No Longer Maintaned HOT 5
- YAML area configuration does not allow to specify the dtype and dtype is ignored in equality comparisons HOT 6
- "resample_nearest" GEO swath to eqc has tiny periodic "glitch" in results HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from pyresample.