In hypertidy/grout: Create Raster Tiles

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-"
)

grout

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))

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))

(t1 <- grout(c(44, 30), blocksize = c(12, 12)))
plot(t1)

(t2 <- grout::grout(c(87, 61), blocksize = c(12, 16)))

plot(t2)

We can generate a table of offset indexes, for use in reading from GDAL (say).

tile_index(t2)

See below for generating wk::rct objects from the tile index.

Or just plot the scheme.

plot(t1)

What for?

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))

## 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"]

Now 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))

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)

TODO

[ ] need more helpers for levels of tiles, even generating all overview levels with these helper extents and dangles
[x] remove need for using sp polygons
[x] set tools for cropping that use the index, not spatial extent (i.e. extent(x, x0, x1, y0, y1))
[x] remove use of sp and raster internally for the data structures, just store the information about the grid/s