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

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Pedro Venâncio-2
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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: 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
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: gdalwarp cutline shift between raster and vector in GDAL 2.3.2

Pedro Venâncio-2
Hi Andre,

Given the magnitude of the shift, I suspected that could be due to the lack of +towgs84 parameters of ESRI 102164.

But I does not understand why it worked correctly in GDAL 2.2.2.

Now I see that

C:\>gdalinfo --version
GDAL 2.3.2, released 2018/09/21

C:\>gdalsrsinfo EPSG:102164
PROJ.4 : +proj=tmerc +lat_0=39.66666666666666 +lon_0=-8.131906111111112 +k=1 +x_0=200000 +y_0=300000 +ellps=intl +units=m +no_defs
OGC WKT :
PROJCS["Lisboa_Hayford_Gauss_IGeoE",
    GEOGCS["GCS_Datum_Lisboa_Hayford",
        DATUM["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"]]

while

pedro@pedro-VirtualBox:~$ gdalinfo --version
GDAL 2.2.2, released 2017/09/15

pedro@pedro-VirtualBox:~$ gdalsrsinfo EPSG:102164
ERROR 1: ERROR - failed to load SRS definition from EPSG:102164

I checked and in both versions, this SRS is listed in proj esri file as

# Lisboa Hayford Gauss IGeoE
<102164> +proj=tmerc +lat_0=39.66666666666666 +lon_0=-8.131906111111112 +k=1.000000 +x_0=200000 +y_0=300000 +ellps=intl +units=m  no_defs <>

and not in the epsg file.

So, is there any explanation for the different behaviour between those versions?

Thank you very much!

Best regards,
Pedro





Andre Joost <[hidden email]> escreveu no dia segunda, 17/12/2018 à(s) 13:17:
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
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: gdalwarp cutline shift between raster and vector in GDAL 2.3.2

Even Rouault-2
On lundi 17 décembre 2018 14:25:02 CET Pedro Venâncio wrote:
> Hi Andre,
>
> Given the magnitude of the shift, I suspected that could be due to the lack
> of +towgs84 parameters of ESRI 102164.

What version of PROJ does your GDAL 2.3.2 build use ? PROJ 5 ?
And your GDAL 2.2.2 build: PROJ 4 ?

If so, that might explain things. The ESRI:102164 translated as a proj string
has no +towgs84, whereas EPSG:20790 has.

With PROJ 4 and how GDAL uses it, due to the absence of a +towgs84 in one of
the source/target CRS, no datum shift is performed, which is what is expected.

With PROJ 5 and how GDAL 2.3.2 uses it, a unwanted datum shift is done (it is
a bit like if the ESRI 102164 had an implicit +towgs84=0,0,0).
I've traced this as https://github.com/OSGeo/gdal/issues/1156
So yes defining -s_srs EPSG:20790 -t_srs EPSG:20790 is the best way to make
sure all CRS are equal to avoid any reprojection attempt.

(the fact that your GDAL 2.2.2 build doesn't find EPSG:102164 is another
issue, probably missing GDAL_DATA, and unrelated, since it doesn't need to be
able to resolve it)

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: gdalwarp cutline shift between raster and vector in GDAL 2.3.2

Pedro Venâncio-2
Hi Even,


What version of PROJ does your GDAL 2.3.2 build use ? PROJ 5 ?
And your GDAL 2.2.2 build: PROJ 4 ?

GDAL 2.3.2 uses PROJ 5.2.0.
GDAL 2.2.3 uses PROJ.4 version 493
GDAL 2.2.2 uses PROJ.4 version 480

 
If so, that might explain things. The ESRI:102164 translated as a proj string
has no +towgs84, whereas EPSG:20790 has.

With PROJ 4 and how GDAL uses it, due to the absence of a +towgs84 in one of
the source/target CRS, no datum shift is performed, which is what is expected.

With PROJ 5 and how GDAL 2.3.2 uses it, a unwanted datum shift is done (it is
a bit like if the ESRI 102164 had an implicit +towgs84=0,0,0).
I've traced this as https://github.com/OSGeo/gdal/issues/1156

I think this is exactly the problem!

C:\>gdalinfo --version
GDAL 2.3.2, released 2018/09/21

C:\>echo "0 0 0" | gdaltransform -s_srs "+proj=longlat +ellps=GRS80 +towgs84=100,100,100" -t_srs "+proj=longlat +ellps=GRS80"
0.000898301199977168 0.000904355202316605 100.001573113725


C:\>gdalinfo --version
GDAL 2.2.3, released 2017/11/20

C:\>echo "0 0 0" | gdaltransform -s_srs "+proj=longlat +ellps=GRS80 +towgs84=100,100,100" -t_srs "+proj=longlat +ellps=GRS80"
0 0 0


pedro@pedro-VirtualBox:~$ gdalinfo --version
GDAL 2.2.2, released 2017/09/15

pedro@pedro-VirtualBox:~$ echo "0 0 0" | gdaltransform -s_srs "+proj=longlat +ellps=GRS80 +towgs84=100,100,100" -t_srs "+proj=longlat +ellps=GRS80"
0 0 0


Do you think this can be fixed soon?

Thank you very much Even!

Best regards,
Pedro



_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev