Converting XY coordinates in the Iberian peninsulate to lat, lon?

Question

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?

snowman2 · Answer

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.

Converting XY coordinates in the Iberian peninsulate to lat, lon?

One Answer

Add your own answers!

Ask a Question