Hello to GDAL community. I have one issue, and I would be very grateful if somebody could give me any kind of help.
I have two rasters of the area around Mount Vesuvius<https://www.google.rs/maps/place/40%C2%B048'45.6%22N+14%C2%B024'51.3%22E/@40.840774,14.2281818,10z/data=!4m5!3m4!1s0x0:0x0!8m2!3d40.81266!4d14.414252>: one spans 20 kilometers, and the other one 50 kilometers (these distances are not so important for the problem). Both rasters are in WGS84 geographic coordinate system, and they overlap between each other. Here is the preview of those two WGS84 rasters<https://www.dropbox.com/s/nza0edj0n4kilfi/vesuvius_radius_20_50KM_overlap.jpg?dl=0> in QGIS. When I use QGIS Raster calculator<https://docs.qgis.org/2.2/en/docs/user_manual/working_with_raster/raster_calculator.html> 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 cell values on the part where they overlap. Here is the preview of difference between the two WGS84 rasters<https://www.dropbox.com/s/c19jj6c3a8if1n2/vesuvius_radius_20-50KM.jpg?dl=0>. 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 Azimuthal equidistant rasters exists<https://www.dropbox.com/s/c4mnpm4szr1l4cq/vesuvius_radius_20-50KM_cubic_aeqd.jpg?dl=0>! I used cubic resampling method, but I get similar results with bilienear method. The reason for this is because once vesuvius_radius_20KM.tif and vesuvius_radius_50KM.tif are reprojected, their grids do not overlap anymore, and their cell sizes slightly differ (79.1578,-79.1578 vs 79.1628,-79.1628). 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-50KM_cubic_aeqd_cellsize100m.jpg?dl=0>. Why? 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_aeqd_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"
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
