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

Etienne Tourigny etourigny.dev at gmail.com
Fri Dec 6 10:44:43 PST 2013 

I have attached a patch to the ticket that fixes the problem. Anyone
interested is welcome to test it.

regards, Etienne


On Fri, Dec 6, 2013 at 1:58 PM, Etienne B. Racine <etiennebr at gmail.com>wrote:

> 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 <etourigny.dev at gmail.com>
>
>> 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/0f588ad1/attachment-0001.html>

More information about the gdal-dev mailing list