Error while transforming/reprojecting geometry using OSGEO/OGR/OSR: "RuntimeError: OGR Error: General Error"

Question

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.

Answer 1

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.