Opened 10 years ago

Closed 10 years ago

#5311 closed defect (fixed)

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

Reported by: etiennebr Owned by: etourigny
Priority: normal Milestone:
Component: Algorithms Version: svn-trunk
Severity: normal Keywords: aggregation, average, mode, resampling
Cc:

Description

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

Attachments (4)

4by4.tif (617 bytes ) - added by etiennebr 10 years ago.
4x4 raster to aggregate
2by2.tif (1.1 KB ) - added by etiennebr 10 years ago.
result of gdal average aggregation
ragg.tif (591 bytes ) - added by etiennebr 10 years ago.
result of R aggregation
patch-5311-1.txt (1.1 KB ) - added by etourigny 10 years ago.
patch #1

Download all attachments as: .zip

Change History (12)

comment:1 by etiennebr, 10 years ago

Summary: Aggregation uses pixels outside the aggregation blocsAggregation probably uses pixels outside the aggregation blocs

comment:2 by etourigny, 10 years ago

Component: defaultAlgorithms
Owner: changed from warmerdam to etourigny
Summary: Aggregation probably uses pixels outside the aggregation blocsAggregation (used in the average resampling algorithm) probably uses pixels outside the aggregation blocs

comment:3 by etourigny, 10 years ago

etiennebr - could you please attach a sample dataset?

by etiennebr, 10 years ago

Attachment: 4by4.tif added

4x4 raster to aggregate

by etiennebr, 10 years ago

Attachment: 2by2.tif added

result of gdal average aggregation

by etiennebr, 10 years ago

Attachment: ragg.tif added

result of R aggregation

by etourigny, 10 years ago

Attachment: patch-5311-1.txt added

patch #1

comment:4 by etourigny, 10 years ago

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

comment:5 by Even Rouault, 10 years ago

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)

comment:6 by etourigny, 10 years ago

Even, thanks for your input.

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

comment:7 by Even Rouault, 10 years ago

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

comment:8 by etourigny, 10 years ago

Resolution: fixed
Status: newclosed

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

Note: See TracTickets for help on using tickets.