[gdal-dev] Why does the difference between two overlapping rasters exist after projection in GDAL?

Didier Bernard deedeebernard at hotmail.com
Thu May 12 04:36:46 PDT 2016


Dear Even,

Thank you!!
I haven't tried with the nearest neighbor resampling method, but your advice on setting the error threshold to 0 (-et 0) worked perfectly with cubic resampling method!!! Now there is no difference between the two Azimuthal equidistant projected rasters! Thank you again.


Not related to this solved issue: you already told me that you are not a C# GDAL developer, but I was wondering, do you happen to know if defining the cell size (resolution) and extent is possible through C# GDAL bindings? Basically using the "-tr" and "-te" parameters from Gdal command prompt syntax?

Regards,
Didier
________________________________________
From: Even Rouault <even.rouault at spatialys.com>
Sent: Wednesday, May 11, 2016 10:42 PM
To: gdal-dev at lists.osgeo.org
Cc: Didier Bernard
Subject: Re: [gdal-dev] Why does the difference between two overlapping rasters exist after projection in GDAL?

Didier,
>
> So if I reproject the initial vesuvius_radius_20KM.tif and
> vesuvius_radius_50KM.tif rasters, but by fixing the cell size (to 100
> meters for example) and by fixing the extents, two rasters would now need
> to have exactly the same values in their cells on the central part (where
> they overlap). But they do not!! There is still a difference between them.
> Smaller difference but it still exists. Here is a preview of that
> difference<https://www.dropbox.com/s/gxdbd1l34jo12uk/vesuvius_radius_20-50
> KM_cubic_aeqd_cellsize100m.jpg?dl=0>.
>
> Why?

Perhaps try adding -et 0 to your gdalwarp command line so exact reprojection
is done. By default there's a small tolerance margin of 0.125 pixel that can
speed up significantly, and might perhaps (I'm not sure) explain the small
difference you see here.

Also did you try with nerarest neighbour to check ?

> That's my question and the point of my email.
> Their grids now match each other, and their cellsizes are exactly 100,100
> meters. So why does the difference exist at the part where they overlap?
>
> I am only interested in projecting the WGS84 rasters to Azimuthal
> equidistant, but just for the sake of checking, I tried using Transverse
> Mercator projection (instead of Azimuthal equidistant), and again the
> difference exists.
>
> I am using GDAL 2.0 and Raster Calculator from QGIS 2.8.7 (to calculate the
> difference between rasters).
>
> Is this some sort of bug related with GDAL 2.0?
>
> I would be very grateful if I could get any kind of reply on this issue.
> Thank you in advance.
>
> P.S
>
>
>
> Additional information:
>
> Here are the initial raster files (in WGS84):
> https://www.dropbox.com/s/pudw8vbgil9ti69/vesuvius_radius_20KM.tif?dl=0
> https://www.dropbox.com/s/k30x4c5442g8028/vesuvius_radius_50KM.tif?dl=0
>
> Here are the Azimuthal equidistant projected rasters with fixed cell sizes
> and extents:
> https://www.dropbox.com/s/cmw1uk3qe09k2u2/vesuvius_radius_20KM_cubic_aeqd_
> cellsize100m.tif?dl=0
> https://www.dropbox.com/s/9yiwc2r6rfd6bgd/vesuvius_radius_50KM_cubic_aeqd_
> cellsize100m.tif?dl=0
>
> And the difference between the two upper Azimuthal equidistant projected
> rasters:
> https://www.dropbox.com/s/jprwjpf44rotnva/vesuvius_radius_20-50KM_cubic_ae
> qd_cellsize100m.tif?dl=0
>
>
>
> To project the initial WGS84 raster files to Transverse Mercator
> projection:
>
> gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32633 -r cubic -of GTiff
> "C:/vesuvius_radius_20KM.tif" "C:/vesuvius_radius_20KM_cubic_tm.tif"
>
>
> gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32633 -r cubic -of GTiff
> "C:/vesuvius_radius_50KM.tif" "C:/vesuvius_radius_50KM_cubic_tm.tif"
>
>
> To project the initial WGS84 raster files to Azimuthal equidistant
> projection:
>
> gdalwarp -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_0=40.81266
> +lon_0=14.414252" -r cubic -of GTiff "C:/vesuvius_radius_20KM.tif"
> "C:/vesuvius_radius_20KM_cubic_aeqd.tif"
>
>
> gdalwarp -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_0=40.81266
> +lon_0=14.414252" -r cubic -of GTiff "C:/vesuvius_radius_50KM.tif"
> "C:/vesuvius_radius_50KM_cubic_aeqd.tif"
>
>
> To project the initial WGS84 raster files to Azimuthal equidistant
> projection by fixing cell size and extents:
>
> gdalwarp -overwrite -te -25000 -18000 14000 21000 -tr 100 100 -s_srs
> EPSG:4326 -t_srs "+proj=aeqd +lat_0=40.81266 +lon_0=14.414252 +x_0=0
> +y_0=0 +ellps=WGS84 +units=m +no_defs" -r cubic -of GTiff
> "C:/vesuvius_radius_20KM.tif"
> "C:/vesuvius_radius_20KM_cubic_aeqd_cellsize100m.tif"
>
>
> gdalwarp -overwrite -te -56000 -48000 44000 51000 -tr 100 100 -s_srs
> EPSG:4326 -t_srs "+proj=aeqd +lat_0=40.81266 +lon_0=14.414252 +x_0=0
> +y_0=0 +ellps=WGS84 +units=m +no_defs" -r cubic -of GTiff
> "C:/vesuvius_radius_50KM.tif"
> "C:/vesuvius_radius_50KM_cubic_aeqd_cellsize100m.tif"

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


More information about the gdal-dev mailing list