Abstract tiling schemes.
Given a grid, impose a tiling scheme. The dangle is overlap the tiles impose (in pixels).
We use the term “block” to refer to the size of each tile in pixels.
We use the simplest specification of a raster grid, which is six
numbers: dimension c(ncol, nrow) and extent
c(xmin, xmax, ymin, ymax).
Consider a raster 16 x 12, with 4 x 4 tiling there is no overlap.
library(grout)
grout(c(4, 4), extent = c(0, 16, 0, 12), blocksize = c(4L, 4L))
#> tiles: 1, 1 (x * y = 1)
#> block: 4, 4
#> dangle: 0, 0
#> tile resolution: 16, 12
#> tile extent: 0, 16, 0, 12 (xmin,xmax,ymin,ymax)
#> grain: 16, 12 (4 : x, 4 : y)But, if our raster has an dimension that doesn’t divide neatly into the block size, then there is some tile overlap. This should work for any raster dimension and any arbitrary tile size.
grout(c(15, 13), extent = c(0, 15, 0, 13), blocksize = c(4L, 4L))
#> tiles: 4, 4 (x * y = 16)
#> block: 4, 4
#> dangle: 1, 3
#> tile resolution: 4, 4
#> tile extent: 0, 16, -3, 13 (xmin,xmax,ymin,ymax)
#> grain: 4, 4 (4 : x, 4 : y)(t1 <- grout(c(44, 30), blocksize = c(12, 12)))
#> tiles: 3, 4 (x * y = 12)
#> block: 12, 12
#> dangle: 4, 6
#> tile resolution: 12, 12
#> tile extent: 0, 48, -6, 30 (xmin,xmax,ymin,ymax)
#> grain: 12, 12 (12 : x, 12 : y)
plot(t1)
(t2 <- grout::grout(c(87, 61), blocksize = c(12, 16)))
#> tiles: 4, 8 (x * y = 32)
#> block: 12, 16
#> dangle: 9, 3
#> tile resolution: 12, 16
#> tile extent: 0, 96, -3, 61 (xmin,xmax,ymin,ymax)
#> grain: 12, 16 (12 : x, 16 : y)
plot(t2)
We can generate a table of offset indexes, for use in reading from GDAL (say).
tile_index(t2)
#> # A tibble: 32 × 11
#> tile offset_x offset_y tile_col tile_row ncol nrow xmin xmax ymin ymax
#> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0 1 1 12 16 0 12 45 61
#> 2 2 12 0 2 1 12 16 12 24 45 61
#> 3 3 24 0 3 1 12 16 24 36 45 61
#> 4 4 36 0 4 1 12 16 36 48 45 61
#> 5 5 48 0 5 1 12 16 48 60 45 61
#> 6 6 60 0 6 1 12 16 60 72 45 61
#> 7 7 72 0 7 1 12 16 72 84 45 61
#> 8 8 84 0 8 1 3 16 84 87 45 61
#> 9 9 0 16 1 2 12 16 0 12 29 45
#> 10 10 12 16 2 2 12 16 12 24 29 45
#> # … with 22 more rowsSee below for generating wk::rct objects from the tile index.
Or just plot the scheme.
plot(t1)
This gives us fine control over the exact nature of the data we can read from large sources.
Consider this large image online:
url <- "https://services.ga.gov.au/gis/rest/services/Topographic_Base_Map/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=Topographic_Base_Map,tilematrixset=default028mm"
dsn <- sprintf("WMTS:%s,layer=Topographic_Base_Map,tilematrixset=default028mm", url)
info <- vapour::vapour_raster_info(dsn)
(overviews <- matrix(info$overviews, ncol = 2, byrow = TRUE))
#> [,1] [,2]
#> [1,] 205851559 165517805
#> [2,] 102925779 82758902
#> [3,] 51462890 41379451
#> [4,] 25731445 20689726
#> [5,] 12865722 10344863
#> [6,] 6432861 5172431
#> [7,] 3216431 2586216
#> [8,] 1608215 1293108
#> [9,] 804108 646554
#> [10,] 402054 323277
#> [11,] 201027 161638
#> [12,] 100513 80819
#> [13,] 50257 40410
#> [14,] 25128 20205
#> [15,] 12564 10102
#> [16,] 6282 5051
#> [17,] 3141 2526
#> [18,] 1571 1263
#> [19,] 785 631
#> [20,] 393 316
#> [21,] 196 158
## level 16 is a reasonable size to experiment with
dm <- overviews[16L, , drop = TRUE]The raster readers in terra or stars gives us an object that can be operated with, if we crop it only those cell values are read - but we have no idea about the underlying tiling of the data source itself.
With GDAL more directly we can find the underlying tile structure, which tells us about the blocked tiling scheme.
info["block"]
#> $block
#> [1] 128 128Now with grout we can actually generate the tile scheme and work with it, let’s say we know that we want a region near tile number 1583. Using the information about the tiles we can find the adjacent tile cell numbers, then use that to crop the original source.
(tile0 <- grout(dm, info$extent, blocksize = info$block))
#> tiles: 40, 50 (x * y = 2000)
#> block: 128, 128
#> dangle: 118, 69
#> tile resolution: 156545.3, 156549.3
#> tile extent: 11148027, 18975292, -6527909, -265936.2 (xmin,xmax,ymin,ymax)
#> grain: 156545.3, 156549.3 (128 : x, 128 : y)
index <- tile_index(tile0)
polys <- wk::rct(index$xmin, index$ymin, index$xmax, index$ymax)
cl <- 1584
tile_dim <- c(max(index$tile_col), max(index$tile_row))
row <- vaster::row_from_cell(tile_dim, cl)
col <- vaster::col_from_cell(tile_dim, cl)
rowcol <- expand.grid(c(-1, 0, 1) + row,
c(-1, 0, 1) + col)
## adjacent (this really needs some work between grout and vaster to make it obvious and easy)
cells <- vaster::cell_from_row_col(tile_dim, rowcol[,1], rowcol[,2])
#cells <- raster::adjacent(tile0$tileraster, 6500, include = TRUE, directions = 8)[, "to"]
exs <- index[cells, c("xmin", "xmax", "ymin", "ymax")]
ex <- c(min(exs$xmin), max(exs$xmax), min(exs$ymin), max(exs$ymax))
## *3 because we got adjacent tiles by row,col above
im <- whatarelief::imagery(source = dsn, extent = ex, dimension = info$block * 3, projection = info$projection)
plot(polys, col = seq_along(polys) == cl)
ximage::ximage(im, extent = ex, add = TRUE)
Finally, we have an in-memory raster of the original source data for very specific tiles.
ximage::ximage(im, extent = ex, asp = 1)
plot(polys[cells], add = TRUE)
text(wk::wk_coords(geos::geos_centroid(polys[cells ]))[, c("x", "y")], lab = cells, col = "firebrick", cex = 1.5)
- need more helpers for levels of tiles, even generating all overview levels with these helper extents and dangles
- remove need for using sp polygons
- set tools for cropping that use the index, not spatial extent (i.e. extent(x, x0, x1, y0, y1))
- remove use of sp and raster internally for the data structures, just store the information about the grid/s
Please note that the grout project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.