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

[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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20181217/b026be8b/attachment.html>