[gdal-dev] how to aggregate an image to a lower resolution using gdal (gdalwarp?)
Etienne B. Racine
etiennebr at gmail.com
Fri Dec 6 07:58:11 PST 2013
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/4d5289ad/attachment-0001.html>
More information about the gdal-dev
mailing list