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

Mon Dec 2 07:41:38 PST 2013

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 <etiennebr at 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.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
>>
>
>
> _______________________________________________
> 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/7197a164/attachment.html>