-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.Rmd
152 lines (100 loc) · 5.06 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-"
)
```
<!-- badges: start -->
[](https://www.tidyverse.org/lifecycle/#experimental)
[](https://CRAN.R-project.org/package=grout)
[](https://github.com/hypertidy/grout/actions)
[](https://github.com/hypertidy/grout/actions/workflows/R-CMD-check.yaml)
<!-- badges: end -->
# 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.
```{r 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.
```{r dangle}
grout(c(15, 13), extent = c(0, 15, 0, 13), blocksize = c(4L, 4L))
```
```{r grid}
(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).
```{r offsets}
tile_index(t2)
```
See below for generating `wk::rct` objects from the tile index.
Or just plot the scheme.
```{r plot}
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:
```{r large}
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.
```{r tile}
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.
```{r scheme}
(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.
```{r specific}
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
---
## Code of Conduct
Please note that the grout project is released with a [Contributor Code of Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms.