[gdal-dev] Difference between ogr2ogr and Python reprojection

Even Rouault even.rouault at spatialys.com
Tue Mar 14 11:25:54 PDT 2023


Nicolas,

ogr2ogr uses the OGRGeometryFactory::transformWithOptions() method that 
has logic to deal with poles, antimerdian, etc. whereas 
OGRGeometry::transform() does basic point-to-point reprojection

You can achieve what you want by using the ogr.GeomTransformer(ct) 
method that returns a ogr.GeomTransformer object that has a 
Transform(ogr.Geometry) method. See below test from GDAL's autotest suite:

def test_ogr_geomtransfomer_default():

     src = osr.SpatialReference()
     src.ImportFromEPSG(32660)

     dst = osr.SpatialReference()
     dst.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
     dst.ImportFromEPSG(4326)

     ct = osr.CoordinateTransformation(src, dst)
     geom = ogr.CreateGeometryFromWkt(
         "LINESTRING(832864.275023695 0,835092.849076364 0)"
     )
     transformer = ogr.GeomTransformer(ct)

     geom_dst = transformer.Transform(geom)
     assert (
         geom_dst.ExportToWkt()
         == "MULTILINESTRING ((179.99 0.0,180 0),(-180 0,-179.99 0.0))"
     )

Even

Le 14/03/2023 à 19:15, Dion-Degodez,Nicolas (ECCC) a écrit :
>
> _Context_
>
> We are trying to reproject a GeoJSON polygon to EPSG:4326.
>
> In some cases, depending on whether we use the ogr2ogr command line or 
> the equivalent in python, the resulting polygon is different.
>
> In those situations, the command line gives us the appropriate 
> polygon, whereas the python code gives us a distorted polygon.
>
> _Question_
>
> Any idea why there is a difference between the command line and the 
> python results? Is there a mistake in our interpretation of the 
> code/how we converted it in python?
>
> _Sample GeoJSON_: https://pastebin.com/hUiVVjFv
>
> Proj4: "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=264 +x_0=0 +y_0=0 
> +R=6371229 +units=m +no_defs"
>
> _Command line that works_: ogr2ogr -s_srs "+proj=stere +lat_0=90 
> +lat_ts=60 +lon_0=264 +x_0=0 +y_0=0 +R=6371229 +units=m +no_defs" 
> -t_srs EPSG:4326 result_command_line.json geojson.json
>
> Result: https://pasteboard.co/cTpZBRamEYv8.jpg
>
> _Python code_ (we would expect the same output as the command line): 
> https://pastebin.com/QgWDBL1W
>
> Result: https://pasteboard.co/8iSv0vbnQlHw.jpg
>
> Thank you!
>
> Nicolas D.
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230314/826099ff/attachment.htm>


More information about the gdal-dev mailing list