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 <[email protected]> > 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 <[email protected]>wrote: > >> 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
