Jan, can you please add your sample dataset and comments in that ticket? Etienne
On Tue, Dec 3, 2013 at 6:37 PM, Etienne B. Racine <etienn...@gmail.com>wrote: > I added a ticket http://trac.osgeo.org/gdal/ticket/5311 > > Etienne > > > 2013/12/2 Etienne Tourigny <etourigny....@gmail.com> > >> The "average" resampling mode of gdalwarp does "average resampling, >> computes the average of all non-NODATA contributing pixels". >> >> It was meant to compute the average of all the pixels in the "aggregation >> window". However, it may have issues in the corners. >> >> I am the author of the average and mode algorithms, and I basically >> copied the aggregation logic from the other algorithms (i.e. which pixels >> are selected for the aggregation), so there may be something wrong in the >> logic. >> >> Certainly, looking at this simple example shows that something is wrong. >> >> I tested average and mode resampling with a fairly large dataset and did >> not find these problems. >> >> Can you create a new bug in the tracker and upload the scripts and data >> used? I don't have much (any) time to work on this but would be happy to >> apply a working patch. >> >> Cheers, Etienne >> >> >> On Mon, Dec 2, 2013 at 12:22 PM, Etienne B. Racine >> <etienn...@gmail.com>wrote: >> >>> I've tried to dig this a bit. I couldn't understand the logic behind >>> gdal aggregation (or downsampling). I've simplified your example using a >>> smaller raster and deterministic values. Maybe someone could enlighten us >>> by looking at the aggregation values ? Note that the lower right cells >>> values were identical in all dimensions I've tried (about 4 - 10). >>> >>> The example is run on GDAL 1.11dev, released 2013/04/13 >>> >>> Thanks, >>> Etienne >>> >>> # in R: >>> >>> require(raster) >>> filei <- '10by10.tif' >>> fileo <- '5by5.tif' >>> >>> dm = 4 >>> r <- raster(matrix(1:(dm^2), dm, dm)) >>> >>> >>> writeRaster(r, filename=filei, overwrite = TRUE) >>> >>> ## aggregate using R aggregate function using the mean >>> r1 <- aggregate(r, fact=2, fun = mean, na.rm = TRUE) >>> >>> file.remove(fileo) >>> cmd <- paste("gdalwarp -r average -tr", paste(res(r1), collapse = " "), >>> filei, fileo) >>> system(cmd) >>> r2 <- raster(fileo) >>> >>> >as.matrix(r) >>> [,1] [,2] [,3] [,4] >>> [1,] 1 5 9 13 >>> [2,] 2 6 10 14 >>> [3,] 3 7 11 15 >>> [4,] 4 8 12 16 >>> > lapply(list(r1, r2, r2-r1), as.matrix) >>> [[1]] >>> [,1] [,2] >>> [1,] 3.5 11.5 >>> [2,] 5.5 13.5 >>> >>> [[2]] >>> [,1] [,2] >>> [1,] 6.0 12.0 >>> [2,] 7.5 13.5 >>> >>> [[3]] >>> [,1] [,2] >>> [1,] 2.5 0.5 >>> [2,] 2.0 0.0 >>> >>> >>> 2013/8/27 Verbesselt, Jan <jan.verbess...@wur.nl> >>> >>> Hi all, >>>> >>>> I have been testing gdalwarp to aggregate (using -r average) an image. >>>> In order to better understand what is happening I have created a >>>> reproducible example within an R environment and compared it with the >>>> aggregate function of the R raster package (see below). There are some >>>> differences between the gdalwarp raster (r2) and the aggregated raster >>>> (r1). >>>> >>>> How is the gdalwarp -r average working? Which pixels are selected for >>>> averaging (e.g.the corner pixels, center pixels, or all within the >>>> aggregation window)? >>>> >>>> If there is a gdal aggregate option to average pixels comparable to the >>>> aggregate function in the R raster package, it would be great as that >>>> would potentially faster, and you could also reproject and aggregate at >>>> once. >>>> >>>> Thanks! >>>> Jan >>>> >>>> http://bfast.r-forge.r-project.org >>>> http://goo.gl/1mBC5F >>>> >>>> >>>> ## R script >>>> require(raster) >>>> filei <- '10by10.tif' >>>> fileo <- '5by5.tif' >>>> >>>> set.seed(123) >>>> r <- raster(ncols=36, nrows=18) >>>> r <- setValues(r, round(runif(ncell(r))*10)) >>>> r >>>> plot(r) >>>> writeRaster(r, filename=filei, overwrite = TRUE) >>>> >>>> ## aggregate using R aggregate function using the mean >>>> r1 <- aggregate(r, fact=2, fun = mean, na.rm = TRUE) >>>> >>>> cmd <- paste("gdalwarp -r average -tr 20 20 -te -180 -90 180 90 ", >>>> filei, " ", fileo, sep = "") >>>> system(cmd) >>>> r2 <- raster(fileo) >>>> >>>> ## compare >>>> plot(r1) >>>> plot(r2) >>>> r1 >>>> r2 >>>> getValues(r1) >>>> getValues(r2) >>>> >>>> ## >>>> plot(r1-r2) >>>> sessionInfo() >>>> >>>> R> sessionInfo() >>>> R version 3.0.1 (2013-05-16) >>>> Platform: x86_64-pc-linux-gnu (64-bit) >>>> >>>> locale: >>>> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=C >>>> LC_COLLATE=C LC_MONETARY=C >>>> [6] LC_MESSAGES=C LC_PAPER=C LC_NAME=C >>>> LC_ADDRESS=C LC_TELEPHONE=C >>>> [11] LC_MEASUREMENT=C LC_IDENTIFICATION=C >>>> >>>> attached base packages: >>>> [1] stats graphics grDevices utils datasets methods base >>>> >>>> other attached packages: >>>> [1] rgdal_0.8-10 raster_2.1-49 sp_1.0-11 >>>> >>>> loaded via a namespace (and not attached): >>>> [1] grid_3.0.1 lattice_0.20-23 tools_3.0.1 >>>> >>>> >>>> rgdal: version: 0.8-10, (SVN revision 478) >>>> Geospatial Data Abstraction Library extensions to R successfully loaded >>>> Loaded GDAL runtime: GDAL 1.10.0, released 2013/04/24 >>>> Path to GDAL shared files: /usr/share/gdal/1.10 >>>> Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480] >>>> Path to PROJ.4 shared files: (autodetected) >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> _______________________________________________ >>>> gdal-dev mailing list >>>> gdal-dev@lists.osgeo.org >>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev >>>> >>> >>> >>> _______________________________________________ >>> gdal-dev mailing list >>> gdal-dev@lists.osgeo.org >>> http://lists.osgeo.org/mailman/listinfo/gdal-dev >>> >> >> >
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev