GRTSmh.Rmd

---
title: "Handling the GRTS master grid"
date: '`r paste("Version",lubridate::now())`'
output:
  html_notebook:
    number_sections: yes
    code_folding: show
    includes:
      in_header: ../header.html
    toc: yes
    toc_float:
      collapsed: no
      smooth_scroll: no
---

```{r setup, message=FALSE, echo=FALSE}
options(stringsAsFactors = FALSE)
library(sp)
library(sf)
library(raster)
library(tidyverse)
library(stringr)
library(n2khab)
# library(plotly)
library(rasterVis)
library(stars)
library(units)
library(tmap)
library(knitr)
opts_chunk$set(
  echo = TRUE,
  dpi = 300
)
```

# Correcting a flaw in the original version


```{r}
GRTSmaster_habitats <- 
    raster("~/Downloads/GRTSmaster_habitats.tif") # i.e. version 1
GRTSmaster_habitats
```

```{r}
inMemory(GRTSmaster_habitats)
```

```{r}
dataType(GRTSmaster_habitats)
```

Are all values unique?

```{r}
GRTSmaster_habitats[] %>% 
    enframe %>% 
    count(value) %>% 
    filter(n > 1)
```

Watch the zeroes!

```{r}
spplot(GRTSmaster_habitats)
```


Oops, see right upper corner. We need to correct this data layer...

First check that the zero value does indeed only occur in the upper right corner:

```{r}
GRTSmaster_habitats[GRTSmaster_habitats == 0,
                    drop = FALSE] %>% 
    spplot
```

Yes, this seems to be a useful criterium.

<!-- ```{r} -->
<!-- flanders <-  st_read("../../n2khab_data/10_raw/flanders") -->
<!-- flanders_500 <-  -->
<!--     flanders %>%  -->
<!--     st_buffer(500) -->
<!-- ``` -->


<!-- ```{r include=FALSE} -->
<!-- # Inspired by code from <https://datacarpentry.org/r-raster-vector-geospatial/08-vector-plot-shapefiles-custom-legend/index.html>, -->
<!-- # however this was too heavy for plotting: -->
<!-- ggplot() + -->
<!--      geom_raster(data = GRTSmaster_habitats %>%  -->
<!--                      as.data.frame(xy = TRUE),  -->
<!--                  aes(x = x,  -->
<!--                      y = y,  -->
<!--                      fill = GRTSmaster_habitats)) + -->
<!--      geom_sf(data = flanders) + -->
<!--      geom_sf(data = flanders_500, -->
<!--              fill = NA, -->
<!--              colour = "red") + -->
<!--     coord_sf() -->
<!-- # compare with: -->
<!-- gplot(GRTSmaster_habitats) + -->
<!--     geom_tile(fill = value) # however, further extension seems difficult here -->
<!-- ``` -->

```{r}
GRTSmaster_simplified <- 
    GRTSmaster_habitats[GRTSmaster_habitats > 0,
                        drop = FALSE]
```

```{r}
GRTSmaster_simplified[] %>% 
    enframe %>% 
    count(value) %>% 
    filter(n > 1)
```



```{r}
GRTSmaster_simplified %>% spplot
```

How many cells have values, and how many have no value?

```{r}
valueslogical <- 
    !(GRTSmaster_simplified %>% 
    values %>% 
    is.na)
valueslogical %>% 
    table
```

Let's compare the 'TRUE' number with the original data to check completeness:

```{r}
valueslogical_orig <- 
    (GRTSmaster_habitats %>% 
    values) > 0
valueslogical_orig %>% 
    sum(na.rm = TRUE)
```

OK, there is a complete match.

Range of the values:

```{r}
GRTSmaster_simplified %>% 
    values %>%
    summary
```


Let's write the corrected version as a GeoTIFF.

```{r}
GRTSmaster_simplified %>% 
    writeRaster(filename = "../../n2khab_data/20_processed/GRTSmaster_habitats.tif",
                datatype = "INT4S")
```

And now let's check the raster header and its summary, from the newly created GeoTIFF:

```{r}
GRTSmaster_habitats2 <- 
    raster("../../n2khab_data/20_processed/GRTSmaster_habitats.tif")
GRTSmaster_habitats2
```


Range of the values:

```{r}
GRTSmaster_habitats2 %>% 
    values %>%
    summary
```

Datatype:

```{r}
dataType(GRTSmaster_habitats2)
```


```{r}
GRTSmaster_habitats2 %>% 
    spplot
```

Are all values unique?

```{r}
GRTSmaster_habitats2[] %>% 
    enframe %>% 
    count(value) %>% 
    filter(n > 1)
```

Do the values match between version 1 and version 2?

```{r}
result <- 
    overlay(x = GRTSmaster_habitats, 
            y = GRTSmaster_habitats2, 
            filename = file.path(tempdir(), "result.tif"),
            datatype = "INT2S",
            fun = function(x, y) {
                return(as.integer(x == y))
            }
    )
```

```{r}
result
```

```{r}
result[] %>% 
    table
```

OK, there is a complete match.

```{r}
spplot(result)
```


# Working on the GRTS master grid

```{r}
GRTSmaster_habitats <- 
    raster("../../n2khab_data/10_raw/GRTSmaster_habitats/GRTSmaster_habitats.tif")
```

## How to make a brick


```{r}
GRTSmh_brick <- 
    brick("../../n2khab_data/10_raw/GRTSmaster_habitats/GRTSmaster_habitats.tif")
```

```{r}
inMemory(GRTSmh_brick)
```


```{r}
GRTSmh_brick2 <- 
    brick(list(GRTSmaster_habitats,
            GRTSmaster_habitats %>% 
                reclassify(rcl = c(1, 1e+9, 1),
                           right = FALSE)
        )
    ) %>% 
    writeRaster(filename = file.path(tempdir(), "GRTSmh_brick2.tif"),
                format = "GTiff", 
                overwrite = TRUE)
names(GRTSmh_brick2) <- (c("decimal", "base4"))
GRTSmh_brick2
```


```{r}
inMemory(GRTSmh_brick2)
```

<!-- ```{r} -->
<!-- writeRaster(GRTSmh_brick2, -->
<!--             filename = file.path(tempdir(),"GRTSmh_brick2.tif"), -->
<!--             datatype = "INT4S", -->
<!--             overwrite = TRUE) -->
<!-- ``` -->

<!-- ```{r} -->
<!-- GRTSmh_brick2 <-  -->
<!--     brick(file.path(tempdir(),"GRTSmh_brick2.tif")) -->
<!-- GRTSmh_brick2 -->
<!-- ``` -->

<!-- ```{r} -->
<!-- names(GRTSmh_brick2)  <-  -->
<!--     c("decimal", "base4") -->
<!-- ``` -->

A way to make an empty brick with `nl` layers, based on a specific layer:

```{r eval=FALSE}
brick(brick(GRTSmaster_habitats), nl = 14)
```

But this is MUCH quicker:

```{r}
brick(extent(GRTSmaster_habitats), 
      nrows = 2843, 
      ncols = 7401, 
      nl = 14, 
      crs = crs(GRTSmaster_habitats))
```



## Going from base 10 to base 4: does it make sense?

We make the function `dec2base4()` based on [baseConvert()](https://github.com/graywh/r-gmisc/blob/master/R/baseConvert.R), which is vectorized and does not give a warning compared to `oro.dicom::dec2base()` (besides, it takes long to load the `oro.dicom` package).
The function gives back the base4 fraction, and probably requires 64 bit operation in R (double is 64 bit; plain integer representations (32-bit) were not sufficient, cf. the non-implemented intermediate `bit64::integer64` base4 representation below).

```{r}
options(scipen = 999,
        digits = 15)
```

```{r }
dec2base4_s <- 
    function(x) {
        ifelse(is.na(x), NA,
            as.double(
                ifelse(x > 0,{
                    d <- floor(log(x, 4) + 1)
                    paste(c("0", "1", "2", "3")[
                        as.integer(abs(diff(x %% 4 ^ seq(d, 0))) %/% 
                                       4 ^ seq(d - 1, 0) + 1)], 
                        collapse = "")
                },
                    '0'
            )) / 10 ^ 13
        )
    }
```


```{r}
dec2base4 <- 
    function(x) sapply(x, dec2base4_s)
```

Examples:

```{r}
dec2base4(NA)
```


```{r}
dec2base4(37784)
```

The maximum base 4 address in the grid:

```{r}
dec2base4(67108857)
```

It works on vectors with NA values:

```{r}
dec2base4(c(14, 15, NA, 456))
```



Let's look at an extracted portion of the GRTSmaster_habitats grid:

```{r}
GRTS_sub <- 
    GRTSmaster_habitats[2000:2050,3000:3050, drop = FALSE]
GRTS_sub %>% 
    spplot
```

```{r}
GRTS_sub_4 <- 
    GRTS_sub %>% 
    calc(function(x){
            dec2base4(x)
        })
GRTS_sub_4 %>% 
    spplot
```

Check whether we can do the inverse and return the original data (again based on [baseConvert()](https://github.com/graywh/r-gmisc/blob/master/R/baseConvert.R)):

```{r}
basetodec <- function(x, base){
    characters <- c(seq(0,9), LETTERS)
    numbers <- structure(seq(0,35), names = characters)    
    x <- strsplit(as.character(x), "")
    l <- lapply(x, length)
    f1 <- function(x, l) {
        sum(numbers[x] * base ^ (seq(l - 1,0)))
    }
    mapply(f1, x, l)
}
```

```{r}
GRTS_sub_4 %>% 
    calc(function(x) basetodec(x * 10 ^ 13, base = 4)) %>% 
    all.equal(GRTS_sub)
```

Olé.

```{r}
stack(GRTS_sub, GRTS_sub_4) %>% 
                    values %>% 
                    as_tibble %>% 
                    as.matrix %>% 
                    head(10)
```

This seems to work perfectly with this raster.

But first, does it make sense spatially; are the original GRTS addresses still reflected in this base4-conversion? Let's extract the 2nd to 6th digit and look at the pattern:

```{r}
max_digits <- 
    c(GRTSmaster_habitats@data@max,
      GRTSmaster_habitats@data@max - 1) %>% 
    dec2base4 %>% 
    as.character %>% 
    (function(x) str_length(x) - 2) %>% 
    max
```

```{r}
extract_digit <- 
    function(x, level, max_digits) {
        withr::with_options(
                  c(scipen = 999), 
                  str_pad(x, max_digits, pad = "0")
                  ) %>% 
            str_sub(level, level) %>% 
            as.numeric
    }
```

```{r}
extract_digit(3333333333320, level = 12, max_digits = max_digits)
```


```{r}
for (i in 2:6) {
GRTS_sub_4 %>% 
    calc(function(x) extract_digit(x * 10 ^ 13, 
                                   level = i, 
                                   max_digits = max_digits)) %>% 
        spplot %>% 
        print
}
```

Yeeha.

We want to be able to aggregate other spatial data at the respective hierarchical levels, and still return this in a spatial way, at the desired level.
Also, at each level we want to have the reverse hierarchical ordering available at that level, so the complete address at that level must be returned.
We will use base 4 **fractions** as address (cf. Stevens & Olsen 2004) as this seems to be more convenient to calculate such an address at each level.

We have two options to implement this:

- do the processing in a dataframe context (group by a specific level, and using simple functions to go to base 4 and calculate the desired grouping level), and then go back to a raster.
However, here we still want to be able not only to use the 32m² raster but also to work with a coarser raster (which is lighter to handle).
- or we do raster calculations only. A brick could even be made that holds the respective levels, but this will still be at the 32 m² level, so quite a burden to calculate.
But that would need to be done only once.

An advantage in both cases is that the coupling between base 10 rank and base 4 addresses is maintained for each raster cell.

So, which is more convenient / feasible?

Update: git-history of this notebook has shown that calculation and storage of such a dataframe is not as efficient as raster calculation.
Those experiments are omitted from the notebook, so we consider raster calculations only.

### Doing raster calculations only

Tried things in a variety of ways, see above and before.
Settling now for the approach that works best.

```{r}
sub10 <- GRTSmaster_habitats[1000:1030, 660:690, drop = FALSE]
sub10 %>% spplot
```

Note what happens with the default datatype 'FLT4S':

```{r}
sub4 <-
    sub10 %>% 
    calc(dec2base4,
        filename = file.path(tempdir(), "sub4.tif"),
        overwrite = TRUE)
```

```{r}
dataType(sub4)
```


```{r}
stack(sub10, sub4) %>% 
                    values %>% 
                    as_tibble %>% 
                    filter(!is.na(layer)) %>% 
                    as.matrix %>% 
                    head(10)
```

Whoopsie, look what happened to the digits!

However when we set the datatype to 'FLT8S', this is solved:

```{r}
sub4 <-
    sub10 %>% 
    calc(dec2base4,
        filename = file.path(tempdir(), "sub4.tif"),
        datatype = "FLT8S",
        overwrite = TRUE)
```


```{r}
sub4 <- raster(file.path(tempdir(), "sub4.tif"))
sub4 %>% spplot
```

```{r}
sub4 %>% dataType()
```

```{r}
sub4[!is.na(sub4)] %>% 
    head(10)
```

```{r}
stack(sub10, sub4) %>% 
                    values %>% 
                    as_tibble %>% 
                    filter(!is.na(layer)) %>% 
                    as.matrix %>% 
                    head(10)
```

This looks perfect and promising!

Let's try to do this for the whole dataset.

```{r}
system.time(
GRTSmh_base4 <-
    GRTSmaster_habitats %>% 
    calc(dec2base4,
        filename = "../../n2khab_data/20_processed/GRTSmh_base4.tif",
        datatype = "FLT8S",
        overwrite = TRUE)
)
```

This took about 15 minutes and the resulting file is 131 MB large, which is rather acceptable, though large.
So we may opt to generate derived layers (cf. evaluation grid) only using decimal ranking numbers, not with the base4 notation.

The file's contents seem fine:

```{r}
GRTSmh_base4 <- raster("../../n2khab_data/20_processed/GRTSmh_base4.tif")
GRTSmh_base4
((GRTSmh_base4 %>% 
    values) > 0) %>% 
    sum(na.rm = TRUE)
GRTSmh_base4[15e+6:(15e+6 + 20)]
```


## Climbing up the GRTS hierarchical level

It feels useful to calculate coarser grids up to the 16 km level, cf.:

```{r}
res <- 32 * 2^(0:12) / 1000
names(res) <- 0:12
res
```

We will declare the names of the above vector (0:12) as the hierarchical levels of the grid.
So we focus on levels 0 to 9.
We want to save a rasterbrick (multilayered tif) with the decimal ranking numbers for those levels.
Saving that as an intermediate dataset will make it easier to work further (dissolving, joining etc).

Let's try this with the subset.

Make an empty brick with the appropriate nlayers, extent, resolution and crs:

```{r warning=FALSE}
sub10brick <- 
    brick(extent(sub4), 
          nrows = 31, 
          ncols = 31, 
          nl = 10, 
          crs = crs(sub4)
          ) %>% 
    writeRaster(filename = file.path(tempdir(), "sub10brick.tif"),
                format = "GTiff", 
                datatype = "INT4S",
                overwrite = TRUE)
names(sub10brick) <- str_c("level", 0:9)
sub10brick
```

Define necessary functions (scalar resp. vectorized) to compute decimal ranks from truncated base-4 fractions:

```{r warning=FALSE}
base4frac2dec_s <- function(x, level) {
            multipliers <- as.matrix(4 ^ ((13 - level - 1):0))
            ifelse(is.na(x), NA, {
                a <- x * 10 ^ level
                a <- round(a - floor(a), 13 - level)
                a <- a %>% 
                    as.character %>% 
                    str_sub(start = 3) %>% 
                    str_pad(width = 13 - level,
                            side = "right",
                            pad = "0") %>% 
                    str_split("", simplify = TRUE) %>% 
                    as.numeric
                t(a) %*% multipliers
            }
            )
}
base4frac2dec <- 
    function(x, level) {
        sapply(x, 
               function(x) base4frac2dec_s(x, level = level))
        }

```

Computation of the brick:

```{r warning=FALSE}
system.time({
for (i in 0:9) {
    sub10brick[[i + 1]] <-
        sub4 %>% 
        calc(function(x) base4frac2dec(x, level = i))
    sub10brick %>%
    writeRaster(filename = file.path(tempdir(), "sub10brick.tif"),
                format = "GTiff",
                datatype = "INT4S",
                overwrite = TRUE)
}
})
```

<!-- ```{r warning=FALSE} -->
<!-- i <- 3 -->
<!-- system.time({ -->
<!-- layer3 <- -->
<!--         sub4 %>%  -->
<!--         calc(function(x) base4frac2dec(x, level = i), -->
<!--              filename = file.path(tempdir(), "layer3.tif"), -->
<!--                 format = "GTiff", -->
<!--                 datatype = "INT4S", -->
<!--                 overwrite = TRUE) -->
<!-- }) -->
<!-- ``` -->


Parallel computation of the brick with `foreach()` may reduce the needed time for calculation, BUT gave different results (inspection of minima and maxima per layer); this may perhaps have to do with a badly exported environment in the foreach function:

```{r warning=FALSE, eval=FALSE}
library(doParallel)
cl <- makeCluster(detectCores())
registerDoParallel(cl)
system.time({
foreach(i = 0:9,
        .export = c("%>%",
                    "calc",
                    "str_sub",
                    "str_pad",
                    "str_split",
                    "writeRaster")) %dopar% {
    sub10brick[[i + 1]] <-
        sub4 %>% 
        calc(function(x) base4frac2dec(x, level = i))
    sub10brick %>%
    writeRaster(filename = file.path(tempdir(), "sub10brick.tif"),
                format = "GTiff",
                datatype = "INT4S",
                overwrite = TRUE)
}
})
# foreach(i = 0:9,
#         .packages = c("raster",
#                       "tidyverse",
#                       "stringr")) %dopar% {
#     sub10brick[[i + 1]] <-
#         sub4 %>% 
#         calc(function(x) base4frac2dec(x, level = i))
#     sub10brick %>%
#     writeRaster(filename = file.path(tempdir(), "sub10brick.tif"),
#                 format = "GTiff",
#                 datatype = "INT4S",
#                 overwrite = TRUE)
# }
# })
stopCluster(cl)
```

Parallel computation of the brick with `clusterR()` ([source](https://stackoverflow.com/questions/44266752/replace-specific-value-in-each-band-of-raster-brick-in-r)) does not work well:

```{r warning=FALSE, eval=FALSE}
beginCluster()
base4frac2dec <- 
    function(x) {
        sapply(x, 
               function(x) base4frac2dec_s(x, level = i))
    }
calcfun <- function(x) calc(x, fun = base4frac2dec)
y <- clusterR(sub4, calc, args = list(fun = base4frac2dec))

# nr_nodes <- 3 #number of cores to use for clusterR function (max recommended: ncores - 1)
# calcfun <- function(x) base4frac2dec(x, level = i)
  #initialize cluster
  beginCluster()
cl <- getCluster()  
library(parallel)
for (i in 0:9) {

  # cat(paste("Currently processing layer:", i + 1,"/",nlayers(sub10brick), "\n"))

  

  # compute
  sub10brick[[i + 1]] <- clusterR(sub4,
                                calc,
                                args = list(fun = base4frac2dec))
    # clusterExport(cl, "i")
    
    # sub10brick[[i + 1]] <- 
    #     sub4 %>% 
    #     calc(function(x){ parSapply(cl, x, base4frac2dec)})
    
}
  #end cluster
  endCluster()
```


Level renaming and control:

```{r warning=FALSE}
names(sub10brick) <- str_c("level", 0:9)
sub10brick
```

Do the same starting from the file:

```{r}
sub10brick <- 
    brick(file.path(tempdir(), "sub10brick.tif"))
names(sub10brick) <- str_c("level", 0:9)
sub10brick
dataType(sub10brick)
```

Having a look at the maps:

```{r}
lapply(as.list(sub10brick), spplot)
```

Great!

# Exploring the created GRTSmh_brick data source

I.e. the datasource generated by code in the `src/generate_GRTS_20_GRTSmh_brick` subfolder.

```{r}
GRTSmh_brick <- 
    brick("../../n2khab_data/20_processed/GRTSmh_brick/GRTSmh_brick.tif")
GRTSmh_brick
```

```{r}
read_GRTSmh_base4frac("../../n2khab_data")
```

```{r}
GRTSmh_brick %>% summary
```




```{r}
read_GRTSmh_base4frac("../../n2khab_data") %>% 
     values %>% 
    max(na.rm = TRUE)
```


```{r}
GRTSmh_brick[[3]] %>% 
     values %>% 
    min(na.rm = TRUE)
```

```{r}
GRTSmh_brick[[10]] %>% 
     values %>% 
    max(na.rm = TRUE)
```

Remind that the cell with the maximum base 4 fraction at level 0 is _not_ necessarily located in the cell with the maximum base 4 fraction at a _higher_ level.

Here, this is the case for level 9 (i.e. the 10th value in below vector):

```{r}
sapply(0:9, function(i) convert_base4frac_to_dec(0.3333333333321, level = i))
```

```{r warning=FALSE}
subbrick <- GRTSmh_brick[1000:1200, 600:800, drop = FALSE]
```


```{r}
lapply(as.list(subbrick[[1:10]]), spplot)
```

```{r}
lapply(as.list(GRTSmh_brick[[7:10]]), spplot)
```


# Make aggregated rasters from a brick

```{r}
GRTSmh_brick <- read_GRTSmh("../../n2khab_data", brick = TRUE)
```


```{r warning=FALSE}
subbrick <- GRTSmh_brick[1000:1200, 600:800, drop = FALSE]
```

After doing some searches, it seems we will have to go via polygonization and rasterization.

## Polygonization step

A method from raster, however this is not efficient:

```{r}
system.time({
sp_layer <- 
    subbrick$layer.6 %>%
    rasterToPolygons(dissolve = TRUE)
})
sp_layer %>% 
    spplot
```

A method from stars, [said](https://github.com/hypertidy/quadmesh#fast-polygonization-of-individual-cells) to be very efficient (cf. the warning in the docs of `spex::polygonize`):

```{r}
system.time({
pol <- 
    subbrick$layer.6 %>% 
    st_as_stars %>% 
    st_as_sf(as_points = FALSE, merge = TRUE)
})
plot(pol)
```

Wow, this is almost 300 times faster!

```{r}
pol
```

Do occurring values coincide with those in the subbrick?

```{r}
pol %>% 
    st_drop_geometry %>% 
    distinct(layer.6) %>% 
    arrange(layer.6) %>% 
    pull(layer.6) %>% 
    all.equal(subbrick$layer.6 %>% 
                  values %>% 
                  unique %>% 
                  sort)
```


## Rasterization step

We could compare `raster::raster()`,`raster::rasterize()`, `fasterize::fasterize()` and `stars::st_rasterize()`.

For that, we need to create a raster template.

However a further option is: `stars::aggregate()`. Maybe first have a look at this; could this prevent the need from going back & forth via polygons?

Trial with `stars::aggregate()` however shows that this does a conversion to the polygon dimension if we supply the polygon as a template.

```{r}
starslayer <- 
    subbrick$layer.6 %>% 
    st_as_stars
aggregate(starslayer, pol, first)
```

Calculating some characteristics of the polygon sf object, excluding non-complete polygons:

```{r}
res <- 
    pol %>% 
    st_cast("MULTILINESTRING") %>% 
    st_length() %>% 
    max / 4
polinterior_coords <- 
    pol %>% 
    mutate(area = st_area(geometry)) %>% 
    filter(area == max(area)) %>% 
    st_geometry %>% 
    st_coordinates
# orig <- 
#     polinterior_coords %>%
#     .[,1:2] %>% 
#     apply(2, min)
ext <- c(xmin = min(polinterior_coords[,1]),
      xmax = max(polinterior_coords[,1]),
      ymin = min(polinterior_coords[,2]),
      ymax = max(polinterior_coords[,2])
      ) %>% 
    extent
```

Plotting the object (i.e. the coloured cells) and the original layout of `pol`):

```{r}
pol %>% 
    mutate(area = st_area(geometry)) %>% 
    filter(area == max(area)) %>% 
    tm_shape() + 
        tm_polygons("layer.6") +
    tm_shape(pol) + 
        tm_polygons(alpha = 0) +
    tm_layout(inner.margins = c(0.2, 0.2, 0.2, 0.2)) +
    tm_legend(show = FALSE)
```

Creating a raster template with 1 row/column added on all sides:

```{r}
rtempl <- 
    raster(ext = ext,
       resolution = drop_units(res),
       crs = crs(subbrick)) %>% 
    extend(c(1, 1))
rtempl
```


### Attempting `stars::st_rasterize()`

Options for GDALRasterizeLayers are [here](https://gdal.org/gdal__alg_8h.html#adfe5e5d287d6c184aab03acbfa567cb1).

```{r}
pol %>% 
    st_buffer(dist = -0.1) %>% # this is done in order to prevent ALL_TOUCHED=TRUE 
                               # from copying the border value of the neighbouring cell
    st_rasterize(rtempl %>% 
                     st_as_stars, 
                 options = "ALL_TOUCHED=TRUE") %>% # in order to not only let cell 
                                                   # centroids determine the value
    write_stars(file.path(tempdir(), "GRTSlowres.tif"))
```

```{r}
GRTSlowres <- 
    raster(file.path(tempdir(), "GRTSlowres.tif"))
GRTSlowres
```

```{r}
GRTSlowres %>% 
    tm_shape() +
    tm_raster() +
    tm_shape(pol) +
    tm_polygons(alpha = 0) +
    tm_text("layer.6", size = 0.6) +
    # tm_layout(inner.margins = c(0.1, 0.1, 0.1, 0.1)) +
    tm_legend(show = FALSE)
```

```{r warning=FALSE}
GRTSlowres[!is.na(GRTSlowres), drop = FALSE] %>% 
    tm_shape() +
    tm_raster() +
    tm_shape(pol) +
    tm_polygons(alpha = 0) +
    tm_text("layer.6", size = 0.6) +
    # tm_layout(inner.margins = c(0.1, 0.1, 0.1, 0.1)) +
    tm_legend(show = FALSE)
```

So this works OK.

### Attempting `raster::raster()`

However from the help, I expect that we will not be able to pass the polygon's attribute value.

```{r}
pol %>%
    raster(ext = ext,
       resolution = drop_units(res),
       crs = crs(subbrick)) %>% 
    values
```

See?

### Attempting `raster::rasterize()`

```{r}
points6 <- subbrick$layer.6 %>% 
  rasterToPoints(spatial = TRUE)
points6 %>% 
  spplot(cex = 0.01)
```

```{r}
points6 %>% 
  rasterize(rtempl,
            field = "layer.6") %>% 
  .[!is.na(.), drop = FALSE] %>% 
  spplot
```

Works perfect, and from applications at the Flemish level (see further and under `src/generate_GRTS_30_GRTSmh_diffres`), I must conclude that `raster::rasterize()` _is very efficient_, at least for points-to-raster calculations.


## Wrapping up

We can also convert stars objects back into raster objects.

```{r}
pol %>% 
    st_buffer(dist = -0.1) %>% # this is done in order to prevent ALL_TOUCHED=TRUE 
                               # from copying the border value of the neighbouring cell
    st_rasterize(rtempl %>% 
                     st_as_stars, 
                 options = "ALL_TOUCHED=TRUE") %>% 
    as("Raster")
```

Yuppekay.

Let's put the preferred chain into one flow.

```{r warning=FALSE}
for (i in 2:8) {
    pol <- 
        subbrick[[i]] %>% 
        st_as_stars %>% 
        st_as_sf(as_points = FALSE, merge = TRUE)
    res <- 
        32 * 2 ^ (i - 1)
    pol_union <- 
        pol %>% 
        st_union
    pol_edge <- 
        pol_union %>% 
        st_buffer(-res) %>% 
        st_difference(pol_union, .)
    polinterior_coords <- 
        pol %>% 
        filter(st_intersects(x = ., 
                             y = pol_edge, 
                             sparse = FALSE)) %>% 
        mutate(area = st_area(geometry)) %>% 
        filter(area == max(area)) %>% 
        st_geometry %>% 
        st_coordinates
    ext <- c(xmin = min(polinterior_coords[,1]),
          xmax = max(polinterior_coords[,1]),
          ymin = min(polinterior_coords[,2]),
          ymax = max(polinterior_coords[,2])
          ) %>% 
        extent
    rtempl <- 
        raster(ext = ext,
               resolution = res,
               crs = crs(subbrick)) %>% 
        extend(c(1, 1))
    colnames(pol)[1] <- "value"
    pol %>% 
        st_buffer(dist = -2) %>% # this is done in order to prevent ALL_TOUCHED=TRUE 
                                   # from copying the border value of the neighbouring cell
                                   # the original value of 0.1 appeared too  small!
        st_rasterize(rtempl %>% 
                         st_as_stars, 
                     options = "ALL_TOUCHED=TRUE") %>% # in order to not only let cell 
                                                       # centroids determine the value
        as("Raster") %>% 
        .[!is.na(.), drop = FALSE] %>% 
        writeRaster(filename = file.path(tempdir(), 
                                         str_c("GRTSmhsub_level", 
                                               i - 1, 
                                               ".tif")),
                    format = "GTiff",
                    datatype = "INT4S",
                    overwrite = TRUE)
    # also storing the pol objects in separate gpkg files (but they could be
    # combined into one, see further):
    pol %>%
        st_write(file.path(tempdir(), 
                                         str_c("pol_level", 
                                               i - 1, 
                                               ".gpkg")),
                 delete_dsn = TRUE
                 )
}
```



```{r}
for (i in 1:7) {
    myraster <- raster(file.path(tempdir(), 
                     str_c("GRTSmhsub_level", 
                           i, 
                           ".tif")))
    mypol <- st_read(file.path(tempdir(), 
                                         str_c("pol_level", 
                                               i, 
                                               ".gpkg"))) 
    (tm_shape(myraster) +
        tm_raster(palette = get_col_regions(),
                  n = 1000) +
        tm_shape(mypol) +
        tm_polygons(alpha = 0) +
        # tm_text("value", size = 0.6, col = "grey") +
        tm_legend(show = FALSE)) %>% 
        print
}
```

This workflow does not work for levels 8 and 9 because the subset of the raster is too small here.



### Geopackage experimenting?

Experimenting with geopackage format (on a layer of `GRTSmh_brick`) to see whether we can put all rasters into one geopackage:

```{r}
r <- GRTSmh_brick$level9
writeRaster(r, filename = file.path(tempdir(), "GRTSmh_diffres.gpkg"),
                    format = "GPKG",
                    datatype = "FLT4S",
                    overwrite = TRUE)
```

The "INT4S" datatype however gives an error, while the "INT2U" and "FLT4S" datatypes do work.
But "INT2U" is useless for the first bands of the brick, so we could use "FLT4S".

HOWEVER, earlier applications on the full dataset have demonstrated a problem: data _do_ change because of the switch from INT4S (in GDAL-terms: Int32) to FLT4S (in GDAL-terms: Float32).
For example, the maximum number 67108857 is changed into 67108856...
So we should stick with INT4S.

The problems with the Int32 are because of the supported data types for GPKG by GDAL itself; see the command `gdalinfo --format GPKG` to see which data types are supported!

So we removed the remainder of the GPKG experiment as it became obsolete.

### Experimenting with Rasterlite

We turn attention to another SQLite implementation, which supports subdatasets _and_ the Int32 datatype.
While it is supported by GDAL and while being an open standard, it is not well supported by R's reading functions.


```{r}
test <-  
  raster(nrows = 6, ncols = 6,
              vals = c(1:35, 67108857)
  )
```


```{r paged.print=FALSE}
test %>% 
  writeRaster(filename = file.path(tempdir(), "test"),
                    format = "GTiff",
                    datatype = "INT4S",
                    overwrite = TRUE)
```

```{r warning=FALSE}
setwd(tempdir())
system("gdal_translate -of Rasterlite -ot Int32 test.tif RASTERLITE:test.sqlite,table=table1")
system("gdal_translate -of Rasterlite -ot Int32 test.tif RASTERLITE:test.sqlite,table=table2")
```

```{r}
gdal_utils(util = "info", 
           file.path(tempdir(), "test.sqlite"))
```

```{r warning=FALSE}
setwd(tempdir())
system("gdal_translate -of GTiff -ot Int32 RASTERLITE:test.sqlite,table=table1 table1.tif")
system("gdal_translate -of GTiff -ot Int32 RASTERLITE:test.sqlite,table=table2 table2.tif")
```

From testing, it appears that `read_stars()` always returns FLT4S data, so it is always converted, which we don't want.
Also, while stars cannot correctly write INT4S, it keeps the numbers in correct shape (`st_as_stars()`), even though it appears that it always converts to FLT4S (I presume there is no difference in memory, for R objects, so it is specific to the raster package to keep track of the file datatype).

We can't read subdatasets from sqlite databases with either `raster()` or `read_stars()`, so we are left with the option of either just storing separate tifs (but note that multipaged tifs are supported from GDAL 3.0.0), or temporarily translating to/from Rasterlite with `gdal_translate` as above.

```{r paged.print=FALSE}
raster(file.path(tempdir(), "table2.tif")) %>% 
  all.equal(test)
```

```{r}
raster(file.path(tempdir(), "table2.tif")) %>% 
  dataType
```

So we could use the following workflow while keeping the Int32 data type:

R raster -> make GeoTIFF -> add into Rasterlite database -> extract subdataset as GeoTIFF -> R raster

### Polygonized versions

As polygonized layers are an intermediate step to the low-resolution raster versions, and because they may come in handy in mapping, it seems valuable to store them as well.

Could this be done using the SQLite driver from GDAL, possibly even using the same SQLite database as above?

```{r}
pol1 <- 
        subbrick[[1]] %>% 
        st_as_stars %>% 
        st_as_sf(as_points = FALSE, merge = TRUE)
pol2 <- 
        subbrick[[2]] %>% 
        st_as_stars %>% 
        st_as_sf(as_points = FALSE, merge = TRUE)
```

```{r}
st_drivers() %>% filter(name == "SQLite")
```

It appears that we cannot use 'SPATIALITE=YES' with layers, so using normal sqlite database:

```{r}
pol1 %>% 
  st_write(file.path(tempdir(), "poltest.sqlite"),
           layer = "pol1",
           driver = "SQLite")
```

```{r}
pol2 %>% 
  st_write(file.path(tempdir(), "poltest.sqlite"),
           layer = "pol2",
           driver = "SQLite")
```

We can read back in:

```{r}
st_read(file.path(tempdir(), "poltest.sqlite"),
           layer = "pol2")
```

Can we add them in the Rasterlite database?

```{r}
pol1 %>% 
  st_write(file.path(tempdir(), "test.sqlite"),
           layer = "pol1",
           driver = "SQLite")
```


```{r}
pol2 %>% 
  st_write(file.path(tempdir(), "test.sqlite"),
           layer = "pol2",
           driver = "SQLite")
```

Yes, we can!
And furthermore, it appears that here the Spatialite standard is used, when inspecting the sqlite database!
For example, the geometry is a binary, as is the case in Spatialite.
It is indicated that all layers have a spatial index, which was missing in the original sqlite db!
Also the original Rasterlite db has several databases referring to Spatialite.

Let's check these layers:

```{r}
st_read(file.path(tempdir(), "test.sqlite"),
           layer = "pol1") %>% 
  all.equal(pol1)
```

A few things seem different, but probably this is only in semantics of metadata:

```{r}
pol1
```

```{r}
pol1t <- st_read(file.path(tempdir(), "test.sqlite"),
           layer = "pol1",
           quiet = TRUE)
pol1t 
```

It seems that the geometry's column name has changed to capitals...

And again check the rasters:

```{r warning=FALSE}
setwd(tempdir())
system("gdal_translate -of GTiff -ot Int32 RASTERLITE:test.sqlite,table=table1 table1.tif")
system("gdal_translate -of GTiff -ot Int32 RASTERLITE:test.sqlite,table=table2 table2.tif")
```

```{r paged.print=FALSE}
raster(file.path(tempdir(), "table2.tif")) %>% 
  all.equal(test)
```

## Upscaling to Flanders: changing my mind


```{r}
dir.create(str_c(datapath, "/20_processed/GRTSmh_diffres"), recursive = TRUE)
```

```{r}
datasetpath <- "../../n2khab_data/20_processed/GRTSmh_diffres"
```


Steps for initialization of a Spatialite database, using the lowest resolution:

```{r warning=FALSE, eval=FALSE}
pol9 <- GRTSmh_brick[[9+1]] %>% 
    st_as_stars %>% 
    st_as_sf(as_points = FALSE, merge = TRUE)
pol7 <- GRTSmh_brick[[7+1]] %>% 
    st_as_stars %>% 
    st_as_sf(as_points = FALSE, merge = TRUE)
res <- 32 * 2 ^ 7
pol_union <-
    pol9 %>%
    st_union
pol_edge <- 
    pol_union %>% 
    st_buffer(-res) %>% 
    st_difference(pol_union, .)
polinterior_coords <- 
    pol7 %>% 
    filter(st_intersects(x = ., 
                         y = pol_edge, 
                         sparse = FALSE)) %>% 
    mutate(area = st_area(geometry)) %>% 
    filter(area == max(area)) %>% 
    st_geometry %>% 
    st_coordinates
ext <- c(xmin = min(polinterior_coords[,1]),
          xmax = max(polinterior_coords[,1]),
          ymin = min(polinterior_coords[,2]),
          ymax = max(polinterior_coords[,2])
          ) %>% 
        extent
rtempl <- 
    raster(ext = ext,
           resolution = res,
           crs = crs(GRTSmh_brick)) %>% 
    extend(c(1, 1))
pol7 %>% 
    st_buffer(dist = -0.1) %>% # this is done in order to prevent ALL_TOUCHED=TRUE 
                               # from copying the border value of the neighbouring cell
    st_rasterize(rtempl %>% 
                     st_as_stars, 
                 options = "ALL_TOUCHED=TRUE") %>% # in order to not only let cell 
                                                   # centroids determine the value
    as("Raster") %>% 
    .[!is.na(.), drop = FALSE] %>% 
    writeRaster(filename = file.path(datasetpath, 
                       "GRTSmh_diffres_level7.tif"),
                format = "GTiff",
                datatype = "INT4S",
                overwrite = TRUE)
setwd(datasetpath)
system("gdal_translate -of Rasterlite -ot Int32 GRTSmh_diffres_level7.tif RASTERLITE:GRTSmh_diffres.sqlite,table=GRTSmh_diffres_level7")
```

Adding the polygon in the sqlite database (also make sqlite separately, for comparison with gpkg):

```{r}
pol7 %>% 
    st_write(file.path(datasetpath, 
                       "GRTSmh_diffres.sqlite"),
           layer = "GRTSmh_polygonized_level7",
           driver = "SQLite")
```

Comparing size to gpkg-version:

```{r}
pol7 %>% 
    st_write(file.path(datasetpath, 
                       "GRTSmh_diffres.gpkg"),
           layer = "GRTSmh_polygonized_level7",
           driver = "GPKG")
```

- The GPKG-version is _much_ smaller than a normal sqlite-database for vector data!
- The TIF file is _much much_ smaller than the Rasterlite version!

So, although it _is possible_ to nicely combine all raster + vector layers into one Rasterlite/Spatialite database file, we will settle with:

- storing the raster layers as separate GeoTIFFs. The format is so much more concise!
- storing the polygonized versions in one GPKG database.

Some reading tests:

```{r}
pol7 <- st_read(file.path(datasetpath, 
                       "GRTSmh_diffres.sqlite"),
           layer = "GRTSmh_polygonized_level7",
           quiet = TRUE)
pol7
```

```{r}
pol7 <- st_read(file.path(datasetpath, 
                       "GRTSmh_diffres.gpkg"),
           layer = "GRTSmh_polygonized_level7",
           quiet = TRUE)
pol7
```

```{r}
setwd(datasetpath)
system("gdal_translate -of GTiff -ot Int32 RASTERLITE:GRTSmh_diffres.sqlite,table=GRTSmh_diffres_level7 GRTSmh_diffres_level7test.tif")
```

```{r}
raster(file.path(datasetpath, "GRTSmh_diffres_level7test.tif"))
```

### Note on the lowest levels (highest resolution)

Using the above workflow (polygonization - rasterization) didn't work out well for level 1 (resolution 64 m), i.e. tests on the result showed repeats of the same value at two different locations.
It's not clear how this became, but it occurs in stars's polygonization step.

Therefore, for the lowests levels we also test a workflow going through points, as the polygonized between-result:

- is more useful as such at higher levels (i.e. for potential applications);
- is not needed at every level for calculating the raster template, because we use the one raster template of level 9 to do rasterization at all levels.

The workflow via points is the `raster::rasterize()` workflow from above.