<div dir="ltr"><div>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).<br>
<br>The example is run on GDAL 1.11dev, released 2013/04/13<br></div><div><br></div><div>Thanks,<br></div>Etienne<br><div><br># in R:<br>require(raster)<br>filei <- '10by10.tif'<br>fileo <- '5by5.tif'<br>
<br>dm = 4<br>r <- raster(matrix(1:(dm^2), dm, dm))<br><br>writeRaster(r, filename=filei, overwrite = TRUE)<br><br>## aggregate using R aggregate function using the mean<br>r1 <- aggregate(r, fact=2, fun = mean, na.rm = TRUE)<br>
<br>file.remove(fileo)<br>cmd <- paste("gdalwarp -r average -tr", paste(res(r1), collapse = " "), filei, fileo)<br>system(cmd)<br>r2 <- raster(fileo)<br><br>>as.matrix(r)<br> [,1] [,2] [,3] [,4]<br>
[1,] 1 5 9 13<br>[2,] 2 6 10 14<br>[3,] 3 7 11 15<br>[4,] 4 8 12 16<br>> lapply(list(r1, r2, r2-r1), as.matrix)<br>[[1]]<br> [,1] [,2]<br>[1,] 3.5 11.5<br>[2,] 5.5 13.5<br>
<br>[[2]]<br> [,1] [,2]<br>[1,] 6.0 12.0<br>[2,] 7.5 13.5<br><br>[[3]]<br> [,1] [,2]<br>[1,] 2.5 0.5<br>[2,] 2.0 0.0<br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">2013/8/27 Verbesselt, Jan <span dir="ltr"><<a href="mailto:jan.verbesselt@wur.nl" target="_blank">jan.verbesselt@wur.nl</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi all,<br>
<br>
I have been testing gdalwarp to aggregate (using -r average) an image.<br>
In order to better understand what is happening I have created a<br>
reproducible example within an R environment and compared it with the<br>
aggregate function of the R raster package (see below). There are some<br>
differences between the gdalwarp raster (r2) and the aggregated raster<br>
(r1).<br>
<br>
How is the gdalwarp -r average working? Which pixels are selected for<br>
averaging (e.g.the corner pixels, center pixels, or all within the<br>
aggregation window)?<br>
<br>
If there is a gdal aggregate option to average pixels comparable to the<br>
aggregate function in the R raster package, it would be great as that<br>
would potentially faster, and you could also reproject and aggregate at<br>
once.<br>
<br>
Thanks!<br>
Jan<br>
<br>
<a href="http://bfast.r-forge.r-project.org" target="_blank">http://bfast.r-forge.r-project.org</a><br>
<a href="http://goo.gl/1mBC5F" target="_blank">http://goo.gl/1mBC5F</a><br>
<br>
<br>
## R script<br>
require(raster)<br>
filei <- '10by10.tif'<br>
fileo <- '5by5.tif'<br>
<br>
set.seed(123)<br>
r <- raster(ncols=36, nrows=18)<br>
r <- setValues(r, round(runif(ncell(r))*10))<br>
r<br>
plot(r)<br>
writeRaster(r, filename=filei, overwrite = TRUE)<br>
<br>
## aggregate using R aggregate function using the mean<br>
r1 <- aggregate(r, fact=2, fun = mean, na.rm = TRUE)<br>
<br>
cmd <- paste("gdalwarp -r average -tr 20 20 -te -180 -90 180 90 ",<br>
filei, " ", fileo, sep = "")<br>
system(cmd)<br>
r2 <- raster(fileo)<br>
<br>
## compare<br>
plot(r1)<br>
plot(r2)<br>
r1<br>
r2<br>
getValues(r1)<br>
getValues(r2)<br>
<br>
##<br>
plot(r1-r2)<br>
sessionInfo()<br>
<br>
R> sessionInfo()<br>
R version 3.0.1 (2013-05-16)<br>
Platform: x86_64-pc-linux-gnu (64-bit)<br>
<br>
locale:<br>
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=C<br>
LC_COLLATE=C LC_MONETARY=C<br>
[6] LC_MESSAGES=C LC_PAPER=C LC_NAME=C<br>
LC_ADDRESS=C LC_TELEPHONE=C<br>
[11] LC_MEASUREMENT=C LC_IDENTIFICATION=C<br>
<br>
attached base packages:<br>
[1] stats graphics grDevices utils datasets methods base<br>
<br>
other attached packages:<br>
[1] rgdal_0.8-10 raster_2.1-49 sp_1.0-11<br>
<br>
loaded via a namespace (and not attached):<br>
[1] grid_3.0.1 lattice_0.20-23 tools_3.0.1<br>
<br>
<br>
rgdal: version: 0.8-10, (SVN revision 478)<br>
Geospatial Data Abstraction Library extensions to R successfully loaded<br>
Loaded GDAL runtime: GDAL 1.10.0, released 2013/04/24<br>
Path to GDAL shared files: /usr/share/gdal/1.10<br>
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]<br>
Path to PROJ.4 shared files: (autodetected)<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
</blockquote></div><br></div>