rspatial / terra

R package for spatial data handling https://rspatial.github.io/terra/reference/terra-package.html
GNU General Public License v3.0
536 stars 89 forks source link

rounding issue in getTileExtents #1564

Closed mdsumner closed 2 months ago

mdsumner commented 2 months ago

There are underlaps and overlaps when getTileExtent(x, SpatRaster) is used, presumably from rounding edge cases. In this example, attempting to break 95x90 into 2x2 tiles results in an overlapping column on the original data, and the same input reduce x4 to 24x23 results in an under-lap of one row.

(I appreciate that breaking into tiles of this nature is somewhat undefined in terms of how alignment and dangle should be handled and it's different for y = SpatRaster than y = numeric, but internal overlaps and underlaps should not occur afaict).

## overplot matrix of extents (cbind(xmin, xmax, ymin, ymax) as from getTileExtents)
rect_extent <- function(x, ...) {
  rect(x[,1], x[,3], x[,2], x[,4], ...)
}
library(terra)
#> terra 1.7.78
f <- system.file("ex/elev.tif", package="terra")
r95x90 <- terra::rast(f)

r2x2 <- rast(r95x90)
dim(r2x2) <- c(2, 2)

e1  <- getTileExtents(r95x90, r2x2)
plot(ext(r95x90), axes = F); rect_extent(e1,  col = rainbow(4, alpha = .5), border = NA) 

e1
#>          xmin     xmax     ymin     ymax
#> [1,] 5.741667 6.141667 49.81667 50.19167
#> [2,] 6.133333 6.533333 49.81667 50.19167
#> [3,] 5.741667 6.141667 49.44167 49.81667
#> [4,] 6.133333 6.533333 49.44167 49.81667

r24x23 <- aggregate(r95x90, 4)
r2x2 <- rast(r24x23)
dim(r2x2) <- c(2, 2)
e2 <- getTileExtents(r24x23, r2x2)
plot(ext(r24x23), axes = F); rect_extent(e2, col = rainbow(4, alpha = .5), border = NA)

e2
#>          xmin     xmax   ymin     ymax
#> [1,] 5.741667 6.141667 49.825 50.19167
#> [2,] 6.141667 6.541667 49.825 50.19167
#> [3,] 5.741667 6.141667 49.425 49.79167
#> [4,] 6.141667 6.541667 49.425 49.79167

Created on 2024-07-18 with reprex v2.0.2

discovered with @njtierney

rhijmans commented 2 months ago

Thanks, the code now ensures that neighboring tiles have shared boundaries.