I added a ticket http://trac.osgeo.org/gdal/ticket/5311
Etienne 2013/12/2 Etienne Tourigny <[email protected]> > 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 <[email protected]>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 <[email protected]> >> >> 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 >>> [email protected] >>> http://lists.osgeo.org/mailman/listinfo/gdal-dev >>> >> >> >> _______________________________________________ >> gdal-dev mailing list >> [email protected] >> http://lists.osgeo.org/mailman/listinfo/gdal-dev >> > >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
