[gdal-dev] Reprojection shapefile using ogr2ogr from EPSG:28992 to EPSG:4326

[Even Rouault]

Le jeudi 22 janvier 2015 19:43:50, Pierric de Laborie a écrit :
> Hello,
> 
> I do not understand why there seems to be a drift (about 100 meters) of the
> polygons contained in a shapefile after reprojecting it to 4326 using the
> following command :
> 
> ogr2ogr -f 'ESRI Shapefile' -wrapdateline -t_srs EPSG:4326
> destinationfile.shp sourcefile.shp
> 
> The original prj file is the following. It is supposed to describe the
> projection EPSG:28992  ( http://spatialreference.org/ref/epsg/28992/  )
> 
> PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",
> SPHEROID["Bessel_1841",6377397.155,299.1528128]],
> PRIMEM["Greenwich",0.0],
> UNIT["Degree",0.0174532925199433]],
> PROJECTION["Double_Stereographic"],
> PARAMETER["False_Easting",155000.0],
> PARAMETER["False_Northing",463000.0],
> PARAMETER["Central_Meridian",5.38763888888889],
> PARAMETER["Scale_Factor",0.9999079],
> PARAMETER["Latitude_Of_Origin",52.15616055555555],
> UNIT["Meter",1.0]]


This definition is outdated and has not the needed datum shift :

$ gdalsrsinfo EPSG:28992

PROJ.4 : '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel 
+towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m 
+no_defs '

OGC WKT :
PROJCS["Amersfoort / RD New",
    GEOGCS["Amersfoort",
        DATUM["Amersfoort",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725],
            AUTHORITY["EPSG","6289"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4289"]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",52.15616055555555],
    PARAMETER["central_meridian",5.38763888888889],
    PARAMETER["scale_factor",0.9999079],
    PARAMETER["false_easting",155000],
    PARAMETER["false_northing",463000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","28992"]]


> 
> Apparently QGIS could understand the projection file correctly and
> reproject it to WGS84 with no problem. But it seems that for some specific
> projections, GDAL does not reproject as expected.

Shapefile .prj files don't include TOWGS84, so that's the likely explanation. 
Try adding -s_srs  EPSG:28992 to your ogr2ogr command line

> 
>  Is it possible to force it by adding the projection specifications in one
> of the files included in the following folder (linux debian install) :
>  /usr/share/gdal   and how ?
> 
> For example, I could see a file named gdal_datum.csv  containing the datum
> "D_Amersfoort" specified in the prj file. However I could not find in any
> other files the projection GCS_Amersfoort or RD_New as specified in the prj
> file.
> 
> Thanks for your help

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com