5

I tried transforming a point from WGS84 (EPSG:4326) to NAD83 / Texas State Mapping System (EPSG:3081) and got a non-descript RuntimeError error. Here's a reproducible example:

from osgeo import ogr, osr

source_osr = osr.SpatialReference()
source_osr.ImportFromEPSG(4326)

target_osr = osr.SpatialReference()
target_osr.ImportFromEPSG(3081)

transform_osr = osr.CoordinateTransformation(source_osr, target_osr)

pt_coords = (-97.7951947301677, 32.7492671878657)

input_point_ogr = ogr.Geometry(ogr.wkbPoint)
input_point_ogr.AddPoint_2D(*pt_coords)

print(input_point_ogr.ExportToIsoWkt())
# POINT (-97.7951947301677 32.7492671878657)

input_point_ogr.Transform(transform_osr)
print(input_point_ogr.ExportToIsoWkt())

Traceback error message:

Traceback (most recent call last):

  File "C:\Users\diasf\AppData\Local\Temp/ipykernel_32464/573320487.py", line 16, in <module>
    input_point_ogr.Transform(transform_osr)

  File "C:\ProgramData\Anaconda3\envs\myenv\lib\site-packages\osgeo\ogr.py", line 7376, in Transform
    return _ogr.Geometry_Transform(self, *args)

RuntimeError: OGR Error: General Error

What's going on here? Why is this happening?

PS: I was able to find out what was going on before posting my question, but I wanted to leave an explicit answer here in case anyone else runs into the same issue in the future. Also, the post in which I found the solution to the issue (this thread) didn't mention the specific error I had.

1
  • 1
    Order of coordinates was changed in GDAL 3.0, according to this
    – Jose
    Commented Feb 15, 2022 at 15:36

1 Answer 1

7

This is because the OSR library expects that geometries in EPSG:4326 will have coordinates ordered as lat-lon:

the Geographic "WGS 84" definition from EPSG, EPSG:4326, (...) uses latitude as the first axis and longitude as the second axis.

Source: OSGEO Wiki (credit to @user30184 who pointed this out in this answer).

In your example, the coordinates are ordered lon-lat. So if you want to want to make the transformation work, you have two options:

  • Invert the order of the coordinates when you create the geometry so that it's lat-lon instead of lon-lat
  • Let OSR know that you are using the lon-lat order instead of the order it expects (which is lat-lon)

Method 1: Inverting the order of the coordinates

from osgeo import ogr, osr

source_osr = osr.SpatialReference()
source_osr.ImportFromEPSG(4326)

target_osr = osr.SpatialReference()
target_osr.ImportFromEPSG(3081)

transform_osr = osr.CoordinateTransformation(source_osr, target_osr)

pt_coords = (-97.7951947301677, 32.7492671878657)

input_point_ogr = ogr.Geometry(ogr.wkbPoint)
input_point_ogr.AddPoint_2D(*pt_coords[::-1])
# Note that, in the command above, we are inverting the order in which we pass 
# the coordinates of the `pt_coords` object to the `AddPoint_2D` function.

print(input_point_ogr.ExportToIsoWkt())
# POINT (-97.7951947301677 32.7492671878657)

input_point_ogr.Transform(transform_osr)
print(input_point_ogr.ExportToIsoWkt())
# POINT (1206249.23361711 1177190.96183372)

Method 2: Informing OSR about the lon-lat order

from osgeo import ogr, osr

source_osr = osr.SpatialReference()
source_osr.ImportFromEPSG(4326)
try: 
    source_osr.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) 
except: 
    pass
# The command above "tells" OSR that the order will be lon-lat.

target_osr = osr.SpatialReference()
target_osr.ImportFromEPSG(3081)

transform_osr = osr.CoordinateTransformation(source_osr, target_osr)

pt_coords = (-97.7951947301677, 32.7492671878657)

input_point_ogr = ogr.Geometry(ogr.wkbPoint)
input_point_ogr.AddPoint_2D(*pt_coords)

print(input_point_ogr.ExportToIsoWkt())
# POINT (-97.7951947301677 32.7492671878657)

input_point_ogr.Transform(transform_osr)
print(input_point_ogr.ExportToIsoWkt())
# POINT (1206249.23361711 1177190.96183372)

Credit for this second method goes to @user30184 and @Simon Gascoin for discussing it in this thread.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.