[gdal-dev] how to aggregate an image to a lower resolution using gdal (gdalwarp?)

Etienne B. Racine etiennebr at gmail.com
Mon Dec 2 06:22:07 PST 2013


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.verbesselt at 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 at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20131202/3e346dd4/attachment.html>


More information about the gdal-dev mailing list