Wrong coordinate transformation with osr.CoordinateTransformation

Question

I'm running a Python script (Python 3.9.16) in ubuntu 22.04.3 LTS. I'm doing it through Docker, althouh I've doubble-checked in Jupyter Notebook and I'm obtaining the same result. I want to transform a bounding box from lat long coordinates (EPSG 4326) to EPSG 32630, but the resulting coordinates are wrong. I'm not getting any errors or warnings.

from osgeo import ogr, osr, gdal

bbox=[-7.84831, 39.71484, -7.73302, 39.81967]

sourceSR = osr.SpatialReference()
sourceSR.ImportFromEPSG(4326) # WGS84 lat lon
targetSR = osr.SpatialReference()
targetSR.ImportFromEPSG(32630)  # WGS84 / UTM zone 30N
coordTrans = osr.CoordinateTransformation(sourceSR, targetSR)
minX_lon, minY_lat, maxX_lon, maxY_lat = bbox
minX, minY, _ = coordTrans.TransformPoint(minX_lon, minY_lat)
maxX, maxY, _ = coordTrans.TransformPoint(maxX_lon, maxY_lat)
trans_bbox = [minX, minY, maxX, maxY]
print(trans_bbox)

The resulting coordinates are: [5696620.047694968, -1177782.6621069708, 5714320.913186516, -1162605.0876848244] and the ones I would expect are: [84341.6581044693, 4407359.8613803, 94846.54347402562, 4418472.2108658375].

I've verified that the X coordinates correspond to longitude and the Y coordinates correspond to latitude. Despite this, I'm still unable to identify what I might be missing. Any insights or suggestions would be greatly appreciated. Thank you!

Answer 1

I came up with the solution!

Before GDAL 3.0 the axis order was longitude first, latitude later. Starting with GDAL 3.0 CRS created with the "EPSG:4326" or "WGS84" strings use the latitude first, longitude second axis order

With oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) you can set the proper order.

I edited my code:

sourceSR = osr.SpatialReference()
sourceSR.ImportFromEPSG(4326) # WGS84 lat lon
sourceSR.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
targetSR.ImportFromEPSG(32630)  # WGS84 / UTM zone 30N
coordTrans = osr.CoordinateTransformation(sourceSR, targetSR)
minX_lon, minY_lat, maxX_lon, maxY_lat = bbox
minX, minY, _ = coordTrans.TransformPoint(minX_lon, minY_lat)
maxX, maxY, _ = coordTrans.TransformPoint(maxX_lon, maxY_lat)
trans_bbox = [minX, minY, maxX, maxY]
print(trans_bbox)

The resulting coordinates are now as expected. I hope this can help someone else.