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 <[email protected]> wrote: > On Thu, Sep 25, 2014 at 8:10 AM, Kyle Shannon <[email protected]> wrote: >> Stephen, >> >> On Wed, Sep 24, 2014 at 2:21 PM, Stephen Roecker >> <[email protected]> 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 >>> [email protected] >>> 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 _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
