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:

enter image description here

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:

enter image description here

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:

enter image description here

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?


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:



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.


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

enter image description here

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
  • Thank you for the reply @AndreJ. Isn't the calculation of the cell value as weighted mean value of the neighboring cells something that is related with the nearest neighbor (-r near in GDAL) resampling method? I am using cubic resampling method, and I also tried with bilinear and got the similar results. – marco May 9 '16 at 13:25
  • All those methods use neighbouring cell values, so you get different results if the new cells are not identical. – AndreJ May 9 '16 at 13:58
  • Thank you @AndreJ. So the reason why cell sizes do not match, is because 50km raster file has a certain number of cells which 20km raster file does not. Due to this new cells, all cells in the 50km raster file are affected, and therefor they can never be the same as the cells of the 20km raster file? Did I understand that correctly? – marco May 9 '16 at 14:27
  • It is not the number of cells, but the cell size. Think of the raster as a polygon grid, that gets deformed from a square to a rhombus during reprojection. – AndreJ May 9 '16 at 15:40
  • Ok, but the cell size depends on the values of the cells, meaning it depends on the number of cells too? – marco May 9 '16 at 16:29

Your Answer

