I have attached a patch to the ticket that fixes the problem. Anyone interested is welcome to test it.
regards, Etienne On Fri, Dec 6, 2013 at 1:58 PM, Etienne B. Racine <etienn...@gmail.com>wrote: > Etienne, I've added the relevant files. I think it's sufficient to debug > the function, but if you think it might be helpful, I could add a larger > raster and the results of the aggregation. > > Etienne > > > 2013/12/6 Etienne Tourigny <etourigny....@gmail.com> > >> 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