Aggregation (used in the average resampling algorithm) probably uses pixels outside the aggregation blocs

[etiennebr]

It seems like the gdal aggregation computes the aggregation statistics on pixels outside of the aggregation blocs. Note that the lower right cells values were identical in all dimensions I've tried (about 4 - 10).

As a comparison, R seems to be doing it fine.

I hypothesize the issue isn't obvious on images because the values are autocorrelated so the artifact is consistent and thus does not appear the the eye.

There's a discussion at ​http://lists.osgeo.org/pipermail/gdal-dev/2013-December/037644.html

The example is run with GDAL 1.11dev, released 2013/04/13

# in R:
# might need install.packages("raster")
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

[etourigny]

etiennebr - could you please attach a sample dataset?

[etiennebr]

4x4 raster to aggregate

[etiennebr]

result of gdal average aggregation

[etiennebr]

result of R aggregation

[etourigny]

patch #1

[etourigny]

etiennebr - can you please apply the attached patch (patch-5311-1.txt) to gdal trunk and test? It works with this dataset and I also looked at the result briefly with a larger dataset (utmsmall.tif - part of the autotest suite).

It would be great if you could test with a large dataset and various scenarios.

thanks

[Even Rouault]

Etienne T., for reference, it appears I introduced, for the other warping algorithms, the 1e-10 in http://trac.osgeo.org/gdal/changeset/22887 for http://trac.osgeo.org/gdal/ticket/4174. Not implying however that it is also necessarily needed for GWKAverageOrModeThread which seems to be particular in the way it works, since for other methods we transform the coordinates at the middle of the pixel (+0.5 shift) whereas for averageOrMode it is done at the upper-left and lower-right corner of the pixels AFAICS, hence a stronger sensitivity to rounding issues if the warping is done just to change resolution (and not reprojection)

[etourigny]

Even, thanks for your input.

So it's ok to remove it? So you see any cases where it would fail?

[Even Rouault]

Yes that's probably OK to remove. It would probably appropriate to add a test to check the case fixed here.

[etourigny]

fixed in trunk (r26696) and 1.10 (r26697) and added an autotest and modified existing tests (r26695)