Geographic Information Systems Asked by Santi Peñate-Vera on April 9, 2021
I have been looking for a way to convert XY coordinates (from a coordinate system based on ETRS89) to latitude and longitude coordinates.
So far I have this function:
from pyproj import Proj, transform
def utm_to_lon_lat(easting, northing):
"""
ETRS89 to latitude, longitude converter
:param easting: Easting coordinates array
:return: latitude, longitude arrays
"""
n = len(easting)
latitude = np.zeros(n)
longitude = np.zeros(n)
inProj = Proj(init='epsg:3857')
outProj = Proj(init='epsg:4326')
for i in range(n):
if not np.isnan(easting[i] + northing[i]):
latitude[i], longitude[i] = transform(inProj, outProj, easting[i], northing[i])
print(easting[i], northing[i], '->', latitude[i], longitude[i])
return latitude, longitude
However, it is not working. The error is that not all the points fall into what I know is the area of the domain, The Iberian peninsula. My guess is that I got the projections wrong, is no it?
The first step I would recommend would be to find the projection that best matches your area of interest. You can do this with the pyproj.CRS class.
Note: the bounds are in this order -> (min_lon, min_lat, max_lon, max_lat)
>>> from pyproj import CRS
>>> CRS("EPSG:25829")
<Projected CRS: EPSG:25829>
Name: ETRS89 / UTM zone 29N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Europe - 12°W to 6°W and ETRS89 by country
- bounds: (-12.0, 34.91, -6.0, 74.13)
Coordinate Operation:
- name: UTM zone 29N
- method: Transverse Mercator
Datum: European Terrestrial Reference System 1989
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
>>> CRS("EPSG:25829").area_of_use.bounds
(-12.0, 34.91, -6.0, 74.13)
>>> CRS("EPSG:25830").area_of_use.bounds
(-6.0, 35.26, 0.0, 80.53)
>>> CRS("EPSG:25831").area_of_use.bounds
(0.0, 37.0, 6.01, 82.41)
>>> CRS("EPSG:25832").area_of_use.bounds
(6.0, 38.76, 12.0, 83.92)
Based on the point: latitude: 40.4369792, longitude: -11.6973331
The best projection looks to be EPSG:25829
Then, you can perform the transformation from latlon (EPSG:4326) to ETRS89 UTM (EPSG:25829) and back again like so:
>>> trans = Transformer.from_crs("EPSG:4326", "EPSG:25829")
>>> trans.transform(40.4369792,-11.6973331)
(271217.25732365483, 4479753.803288419)
>>> trans = Transformer.from_crs("EPSG:25829", "EPSG:4326")
>>> trans.transform(271217.25732365483, 4479753.803288419)
(40.43697919999999, -11.6973331)
Note: axis order can be tricky. If you prefer lon lat (traditional xy order) I would look into using the always_xy
option on the Transformer.
Answered by snowman2 on April 9, 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