Opened 3 years ago

Closed 3 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 3 years ago.
4x4 raster to aggregate
2by2.tif (1.1 KB) - added by etiennebr 3 years ago.
result of gdal average aggregation
ragg.tif (591 bytes) - added by etiennebr 3 years ago.
result of R aggregation
patch-5311-1.txt (1.1 KB) - added by etourigny 3 years ago.
patch #1

Download all attachments as: .zip

Change History (12)

comment:1 Changed 3 years ago by etiennebr

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

comment:2 Changed 3 years ago by etourigny

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 Changed 3 years ago by etourigny

etiennebr - could you please attach a sample dataset?

Changed 3 years ago by etiennebr

Attachment: 4by4.tif added

4x4 raster to aggregate

Changed 3 years ago by etiennebr

Attachment: 2by2.tif added

result of gdal average aggregation

Changed 3 years ago by etiennebr

Attachment: ragg.tif added

result of R aggregation

Changed 3 years ago by etourigny

Attachment: patch-5311-1.txt added

patch #1

comment:4 Changed 3 years ago by 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

comment:5 Changed 3 years ago by 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)

comment:6 Changed 3 years ago by etourigny

Even, thanks for your input.

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

comment:7 Changed 3 years ago by Even Rouault

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

comment:8 Changed 3 years ago by etourigny

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.