[gdal-dev] Fwd: gdalwarp -r average bug?

[Stephen Roecker]

After doing some more testing with gdal 11.0 from the QGIS installer,
I noticed that the gdalwarp -r average function works as advertisted
when applied to a raster dataset using a projected cooridinate system.
The dataset from my example was using a geographic coordinate system.

Stephen

On Thu, Sep 25, 2014 at 10:13 AM, Kyle Shannon <kyle at pobox.com> wrote:
> On Thu, Sep 25, 2014 at 8:10 AM, Kyle Shannon <kyle at pobox.com> wrote:
>> Stephen,
>>
>> On Wed, Sep 24, 2014 at 2:21 PM, Stephen Roecker
>> <stephen.roecker at gmail.com> wrote:
>>> Out of curiousity the other day I compared the results of gdalwarp (-r
>>> average) against  the raster R package aggregate(fun=mean) function
>>> for aggregating a raster to a coarser resolution. I was suprized how
>>> different the results of gdalwarp were from raster. When zooming in
>>> and manually averaging the overlapping cells, the results of gdal were
>>> off.
>>>
>>> The elevation difference between the raster and gdal aggregated
>>> rasters only averaged 0.66 meters, but had a max of 12 meters. Also
>>> when subtracting the raster and gdal aggregated rasters, the resulting
>>> subtracted layer looked like a hillshade, suggesting the GDAL
>>> aggregated raster was shifted. However the raster and GDAL rasters
>>> overlapped perfectly, so I can only assume the shift occured during
>>> the aggregation process.
>>>
>>> Can someone explain gdal's behavior to me? Why the difference, is this
>>> a bug in gdal? gdalwarp claims it's averaging all the overlapping
>>> cells except the NA. That doesn't seem to be the case. FYI I'm using
>>> GDAL 1.10.1
>>>
>>> See a reproduceable R example below.
>>>
>>> Stephen
>>>
>>> library(gdalUtils)
>>> library(raster)
>>>
>>> src_dataset <- system.file("external/tahoe_lidar_bareearth.tif",
>>> package="gdalUtils")
>>> test <- raster(x=src_dataset)
>>>
>>> writeRaster(test, "test.tif", overwrite=T)
>>> gdal_setInstallation(search_path="C:/ProgramData/QGIS/QGISDufour/bin",
>>> rescan=T, verbose=T)
>>> gdalwarp(srcfile=src_dataset, dstfile="test_gdal.tif", of="GTiff",
>>> r="average", ot="Float32", tr=res(test)*3, overwrite=TRUE,
>>> verbose=TRUE)
>>>
>>> aggregate(x=raster(src_dataset), fact=3, filename="test_raster.tif",
>>> format="GTiff", NAflag=-99999,
>>>   progress="text", overwrite=T)
>>> _______________________________________________
>>> gdal-dev mailing list
>>> gdal-dev at lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>> Without looking too far into it, it seems similar to:
>>
>> http://trac.osgeo.org/gdal/ticket/5311
>>
>> --
>> Kyle
>
> Apologies, it does appear to be different, my mistake.
>
> --
> Kyle