Why does the difference between two overlapping rasters exist?

Question

I have two similar rasters of the Mount Vesuvius: one spans 20 kilometers, and the other one 50 kilometers. Both rasters are in WGS84 geographic coordinate system, and they overlap between each other:

When I use QGIS Raster calculator to calculate the difference between them, a new raster file is created which is totally black - so it shows that there is absolutely no difference in pixel values on the part where they overlap:

However, when I project both rasters to Azimuthal equidistant projection, and again check for the difference with Raster calculator, the difference between these two reprojected rasters exists! Here is how it looks like:

Why is this so?

The difference does not exist between the two rasters when they are in WGS84 geographic coordinate system, but once they are projected to Azimuthal equidistant projection, a significant difference appears?

Why?

I used the cubic resampling method, but I got the similar upper difference with the bilinear one.

Just for the sake of checking, I tried using Transverse Mercator projection (instead of Azimuthal equidistant), and again a difference exists between initial two projected rasters. Why is that so?

---

Additional information: To project the initial WGS84 raster files, I used the following syntax:

gdalwarp -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_ts=40.81266 +lon_0=14.414252" -r cubic -of GTiff "C:/vesuvius_radius_20KM.tif" "C:/vesuvius_radius_20KM_cubic_aeqd.tif"

To project them to Transverse Mercator:

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"

Here are the initial raster files (in WGS84):

https://www.dropbox.com/s/98hcftjsnmrqs8p/vesuvius_radius_20KM.tif?dl=0 https://www.dropbox.com/s/a68921e6tpszt0f/vesuvius_radius_50KM.tif?dl=0

Here are their Azimuthal equidistant projected rasters:

https://www.dropbox.com/s/4mexe0rmieam9ri/vesuvius_radius_20KM_cubic_aeqd.tif?dl=0 https://www.dropbox.com/s/0hcle0b151kbod1/vesuvius_radius_50KM_cubic_aeqd.tif?dl=0

And the difference between two Azimuthal equidistant projected rasters:

https://www.dropbox.com/s/2w39hgb9kcrzb8h/vesuvius_radius_20-50KM_cubic_aeqd.tif?dl=0

Answer 1

If you reproject a raster to a different projection, the cell sizes do not match to the original. For this reason, the new cell value will be calculated with respect to neighbouring cells as a weighted mean value.

In your original files, both cell sizes matched exactly (0.000833333 degrees), and the cell values were identical. But after reprojection, cell size (82.7608 m vs 82.7644 m) and values are not identical anymore.

---

UPDATE

To clarify, I have created vector grids from your raster files, with the same extent and cell size:

The green grid is from the original 50km- and 20km-raster, perfectly overlapping with a cell size of 0.000833333° (3 seconds). Red lines are from the 50km-aeqd grid, and blue ones from the 20km-aeqd-grid.

For each red and blue cell, the values from the touched green cells will be taken to calculate the new cell value. Since the red and blue cell do not align, you get different cell values for both grids.

After reprojection, the rasters are not exactly rectangular anymore. GDAL tries to create a raster size almost fitting to the original. To get a full number of rows and columns from the new extent, it gets different spacing for the 50km and 20km gird. You can force the cell size to a certain value with -tr option in gdalwarp, but the grids may still be shifted against each other. You would have to define the target extents with -te as well to make sure that both rasters align.

So something like this should work:

gdalwarp -overwrite -te -28000 4500000 15000 4540000 -tr 100 100 -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_0=0 +lon_0=14.414252 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" -of GTiff vesuvius_radius_20KM.tif 20km-neu.tif
gdalwarp -overwrite -te -60000 4470000 50000 4570000 -tr 100 100 -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_0=0 +lon_0=14.414252 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" -of GTiff vesuvius_radius_50KM.tif 50km-neu.tif