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

[Even Rouault]

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