Geographic Information Systems Asked by Kianwee Chen on July 12, 2021
I am trying to convert between EPSG 4326 and 4978. I have wrote the following codes:
import pyproj
lon = -74.6604194
lat = 40.3450248
transformer = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:4978", always_xy = True)
print(transformer)
xy = transformer.transform(lon, lat)
print('Projected XY', xy)
lonlat = transformer.transform(xy[0], xy[1], direction='INVERSE')
print('Lon Lat', lonlat)
This is the result it produces. As you can see it converts the lon/lat to the right xy. However, when I try to convert xy back to lon/lat, I always get a 0.0 lat value. I have tried with other xy values and I always get a 0.0 lat. What seems to be wrong with my code?
proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step proj=cart ellps=WGS84
Projected XY (1287775.7028415564, -4694570.895022585)
Lon Lat (-74.6604194, 0.0)
I figured out what is wrong. Posting it here hoping it might be useful for others facing the same mistake I made. I think as 4978 is a ecef, you need to input the altitude(Z) to get a proper transformation.
import pyproj
lon = -74.6604194
lat = 40.3450248
alt = 50
transformer = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:4978", always_xy = True)
print(transformer)
xyz = transformer.transform(lon, lat, alt)
print('Projected XY', xyz)
lonlat = transformer.transform(xyz[0], xyz[1], xyz[2], direction='INVERSE')
print('Lon Lat', lonlat)
The results.
proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step proj=cart ellps=WGS84
Projected XY (1287785.7839034656, -4694607.645412987, 4107291.4290166595)
Lon Lat (-74.6604194, 40.3450248, 50.0)
Answered by Kianwee Chen on July 12, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP