[gdal-dev] gdalwarp cutline shift between raster and vector in GDAL 2.3.2

[Andre Joost]

Hi Pedro,

if you compare gdalsrsinfo epsg:20790 and gdalsrsinfo EPSG:102164 from 
ESRI, you see that ESRI omits the datum shift between Lisbon Hayford and 
WGS84.

So you should use -s_srs EPSG:20790 to get the warping correctly.

HTH,
Andre Joost

Am 17.12.18 um 11:09 schrieb Pedro Venâncio:
> Hi,
>
> I'm getting an error when clipping a raster by a mask layer, in GDAL 2.3.2.
>
> This beaviour does not happen in GDAL 2.2.2.
>
> To be more specific, in a Windows machine with GDAL 2.3.2, the clip is
> done, but the image is clipped about 200 meters away from the correct
> bounding box of the vector layer.
>
> In Linux, with GDAL 2.2.2, the clip is done exactly by the vector bounding
> box, as expected.
>
> The command I'm using is just:
>
> gdalwarp -of GTiff -tr 4.235097023884323 -4.2350970238843235 -tap -cutline
> PATH_TO_FOLDER\Cartograma_20790_19D.gpkg -dstalpha PATH_TO_FOLDER\19-D.jpg
> PATH_TO_FOLDER\19-D_clip.tif
>
> The distance of the error made me think that could be something related to
> the coordinate system. In fact, the aux.xml file from the raster layer (in
> jpg format), and which must have been generated in some ESRI software, has
> the EPSG:20790 coordinate system declared as:
>
>        <SpatialReference xsi:type="typens:ProjectedCoordinateSystem">
>
> <WKT>PROJCS["Lisboa_Hayford_Gauss_IGeoE",GEOGCS["GCS_Datum_Lisboa_Hayford",DATUM["D_Datum_Lisboa_Hayford",SPHEROID["International_1924",6378388.0,297.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",200000.0],PARAMETER["False_Northing",300000.0],PARAMETER["Central_Meridian",-8.131906111111112],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",39.66666666666666],UNIT["Meter",1.0],AUTHORITY["Esri",102164]]</WKT>
>          <XOrigin>-5423400</XOrigin>
>          <YOrigin>-14095000</YOrigin>
>          <XYScale>450251902.2805022</XYScale>
>          <ZOrigin>-100000</ZOrigin>
>          <ZScale>10000</ZScale>
>          <MOrigin>-100000</MOrigin>
>          <MScale>10000</MScale>
>          <XYTolerance>0.001</XYTolerance>
>          <ZTolerance>0.001</ZTolerance>
>          <MTolerance>0.001</MTolerance>
>          <HighPrecision>true</HighPrecision>
>          <WKID>102164</WKID>
>          <LatestWKID>102164</LatestWKID>
>        </SpatialReference>
>
> Simply making a gdal_translate -a_srs EPSG:20790 to the original raster
> layer, produces the following info in the aux.xml:
>
>      <PAMDataset>
>        <SRS>PROJCS["Lisbon (Lisbon) / Portuguese National
> Grid",GEOGCS["Lisbon
> (Lisbon)",DATUM["Lisbon_1937_Lisbon",SPHEROID["International
> 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-304.046,-60.576,103.64,0,0,0,0],AUTHORITY["EPSG","6803"]],PRIMEM["Lisbon",-9.131906111111112,AUTHORITY["EPSG","8902"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4803"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",39.66666666666666],PARAMETER["central_meridian",1],PARAMETER["scale_factor",1],PARAMETER["false_easting",200000],PARAMETER["false_northing",300000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","20790"]]</SRS>
>        <GeoTransform>  1.6169619321196713e+05,  4.2462876699946506e+00,
> -9.3136981756810518e-03,  3.6360765297460271e+05,  4.1234595692999079e-03,
> -4.2324755960434057e+00</GeoTransform>
>      </PAMDataset>
>
> Making the gdalwarp cutline to this new raster layer, with EPSG:20790
> assigned, already gives the correct output.
>
> So, the strange thing is the diferent beaviour between GDAL 2.2.2 and
> 2.3.2. Can this be a bug in GDAL 2.3.2?
>
> Here is a sample dataset to test:
>
> https://cld.pt/dl/download/0fcb8923-ed23-45a0-b518-0e5168d0cc0b/19-D.zip
>
> Thank you very much.
>
> Best regards,
> Pedro Venâncio
>
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>