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

Fri Dec 6 07:18:44 PST 2013

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 <etiennebr at gmail.com>wrote:

> I added a ticket http://trac.osgeo.org/gdal/ticket/5311
>
> Etienne
>
>
> 2013/12/2 Etienne Tourigny <etourigny.dev at 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 <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/20131206/40b14b61/attachment.html>