this example for confidence in the algo

[mdsumner]

Needs some fixing for the output extent

---

title: "Minute Detail Rasterizing"
author: "Michael Sumner"
format: html
editor: visual
editor_options:
chunk_output_type: console

There are multiple ways of rasterizing in R. This is a whirlwind tour with a tiny dataset to make sure
an attempt at a more abstracted approach is in alignment with more standard "spatial-ways" of doing things.

I'm not going to explain or excuse the code below, I just want to record it for now.

polstemp <- tempfile(fileext = "fgb")
pols <- silicate::minimal_mesh
sf::write_sf(pols, polstemp, driver = "FlatGeobuf")

## wrap everything in a function so we have just one parameter
fun <- function(n) {
library(silicate)
library(vaster)
library(controlledburn)

sc <- SC0(minimal_mesh)

e0 <- c(apply(sc_vertex(sc), 2, range))
dxy <- diff(e0)[c(1, 3)] / n

res <- round(max(dxy), 3)

e <- vaster::buffer_extent(e0 + c(-1, 1, -1, 1) * rep(dxy, each = 2), res)
dm <- as.integer(diff(e)[c(1, 3)]/res)
print(dm)

rr <- raster::raster(rast(ext(e), ncols = dm[1], nrows = dm[2]))
tr <- rasterize(vect(polstemp), rast(rr))

fr <- fasterize::fasterize(pols, rr)
plot(tr, col = scales::alpha("black", .5))

#par(xpd = NA)
#plot.new(); plot.window(e[1:2], e[3:4], asp = 1)
#plot_extent(e, asp = 1, add = TRUE)
plot(x, add = T)
abline(v = x_corner(dm, e),
       h = y_corner(dm, e), 
       col = "darkgrey")


cb <- controlledburn:::burn_polygon(pols, extent = e,
                               dimension = dm)

## our index is triplets of start,end,line where the polygon edge was detected - 
## this essentially an rle by scanline of start,end polygon coverage
index <- matrix(unlist(cb, use.names = F), ncol = 4L, byrow = TRUE) + 1 ## plus one because 0-index internally


#cb <- controlledburn:::burn_polygon(minimal_mesh, extent = e, dimension = dm)
#str(cb)

## what does controlleburn know?
index <- matrix(unlist(cb, use.names = F), ncol = 4L, byrow = TRUE) + 1 ## plus one because 0-index internally
x0 <- x_from_col(dm, e, index[,1L])
x1 <- x_from_col(dm, e, index[,2L])
y0 <- y_from_row(dm, e, index[,3L])
y1 <- y_from_row(dm, e, index[,3L])

segments(x0, y0, x1, y1)
if (n < 34) {
points(x0, y0, pch = "+", cex = .7)
points(x1, y1, pch = "x", cex = .7)

points(xyFromCell(fr, seq_len(prod(dm))[!is.na(values(fr))]), cex = 2)
} else {
  points(x0, y0, pch = 19, cex = .25, col = "green")
  points(x1, y1, pch = 19, cex = .25, col = "red")

}

invisible(NULL)
}

The parameter n is the number of pixels we want in the plane. So, first let's do a sensible default, a raster that is 10x10.

What do we see here? The yellow is good old terra::rasterize(), a single pixel raster in which was detected our two polygons

i <- 29
i <- i - 1
par(bg = "grey70")
fun(5)

[mdsumner]

more

---
title: "Minute Detail Rasterizing"
author: "Michael Sumner"
format: html
editor: visual
editor_options: 
  chunk_output_type: console
---

There are multiple ways of rasterizing in R. This is a whirlwind tour with a tiny dataset to make sure
an attempt at a more abstracted approach is in alignment with more standard "spatial-ways" of doing things. 

I'm not going to explain or excuse the code below, I just want to record it for now. 


```{r code}
## Two Multipolygons (start point duplicated on each ring)
## a) 8 points, with nested 6 points as a hole 
## b) 5 points
w <- c(a = "MULTIPOLYGON (((0 0, 0 1, 0.75 1, 1 0.8, 0.5 0.7, 0.8 0.6, 0.7 0, 0 0), 
                       (0.2 0.2, 0.5 0.2, 0.5 0.4, 0.3 0.6, 0.2 0.4, 0.2 0.2)))", 
       b = "MULTIPOLYGON (((0.7 0, 0.8 0.6, 1.1 0.63, 
                      1.23 0.3, 0.7 0)))")

## create a spatial type
sp <- wk::as_wkb(w)

sp <- wk::as_wkb(silicate::inlandwaters$geom)
library(wk)
library(vaster)
library(controlledburn)
library(terra)

## every time we rasterize we're going to expand the domain a bit, this is not important we ## just don't want to hug the polygon extremeties for our pictures (we could though)

mkbb <- function(x, f = 0.1) {
  bb <- as.numeric(wk::wk_bbox(x))
  xlim <- bb[c(1, 3)]
  bxlim <- xlim + c(-1, 1) * diff(xlim) * f
  ylim <- bb[c(2, 4)] 
  bylim <- ylim + c(-1, 1) * diff(ylim) * f
  ## xmin, ymin, xmax, ymin
  c(rbind(bxlim, bylim))
}

mkrct <- function(x, f = .1) {
  b <- mkbb(x, f)
  wk::rct(b[1], b[2], b[3], b[4])
}

fun <- function(x, n = 10, f = .1) {
 ## we always work with this extent, we know it exceeds our data (or we could set f = 0)
 bb <- mkbb(x, f = f)
 rct <- mkrct(x, f = f)
 ex <- bb[c(1, 3, 2, 4)]
 ## what size are the pixels
 xlim <- bb[c(1, 3)]
 ylim <- bb[c(2, 4)]
 ## always use equal res (we can be more GDAL-warper later)
 res <- diff(xlim) / n

 dm <- dim(rast(ext(ex), res = res))[2:1]
 ## oldschool
rr <- raster::raster(rast(ext(ex), ncols = dm[1], nrows = dm[2]))
## newschool
tr <- rasterize(vect(as.character(wk::as_wkt(x))), rast(rr))
v <- values(tr)
v[is.nan(v)] <- NA

## fast school (but approximate)
fr <- fasterize::fasterize(x, rr)



cb <- controlledburn:::burn_polygon(sf::st_sf(g = sf::st_as_sfc(x)), extent = ex,
                               dimension = dm)

suppressWarnings(exact <- exactextractr::exact_extract(rr, sf::st_sf(g = sf::st_as_sfc(x)), include_cell = TRUE))

cell <- do.call(rbind, lapply(exact, function(x) as.matrix(x[,2:3])))
cell <- cell[cell[,2] < 1, , drop = FALSE]
rrr <- rast(rr)

rrr[cell[,1]] <- cell[,2]

## our index is triplets of start,end,line where the polygon edge was detected - 
## this essentially an rle by scanline of start,end polygon coverage
index <- matrix(unlist(cb, use.names = F), ncol = 4L, byrow = TRUE) + 1 ## plus one because 0-index internally
plot.new()
plot.window(xlim, ylim, asp = 1, xaxs = "i", yaxs = "i")


plot(rct, lty = 2, add = T)
#plot(rrr, range = c(0, 1), add = TRUE)
box(col = "white", lwd = 3)

#plot(tr, col = scales::alpha("black", .5), add = TRUE, legend = F)

plot(x, add = T, lwd = 3)
# abline(v = x_corner(dm, ex),
#        h = y_corner(dm, ex), 
#        col = "firebrick")
plot(as.polygons(tr, aggregate = F, na.rm = FALSE), col = NA, border = "black", add = TRUE)

fbf <- scales::alpha("firebrick", 0.4)
plot(as.polygons(tr, aggregate = F), 
     col = fbf, border = fbf, add = TRUE)

## what does controlleburn know?
colnames(index) <- c("start", "end", "row", "id")
x0 <- x_from_col(dm, ex, index[,1L])
x1 <- x_from_col(dm, ex, index[,2L])
y0 <- y_from_row(dm, ex, index[,3L])
y1 <- y_from_row(dm, ex, index[,3L])

segments(x0, y0, x1, y1)
if (n < 34) {
points(x0, y0, pch = "|", cex = .7)
points(x1, y1, pch = ">", cex = .7)

#points(xyFromCell(fr, seq_len(prod(dm))[!is.na(values(fr))]), cex = 2)
} else {
  points(x0, y0, pch = 19, cex = .15, col = "green")
  points(x1, y1, pch = 19, cex = .15, col = "red")

}
title("fasterize", line = -4)

plot.new()
plot.window(xlim, ylim, asp = 1, xaxs = "i", yaxs = "i")

plot(rct, lty = 2, add = T)

box(col = "white", lwd = 3)
plot(as.polygons(tr, aggregate = F, na.rm = FALSE), col = NA, border = "black", add = TRUE)

plot(as.polygons(tr, aggregate = F), col = fbf, border = fbf, add = TRUE)

plot(rrr, col = grey.colors(100, 0, 1, rev = TRUE), range = c(0, 1), add = TRUE)
plot(x, add = TRUE, lwd = 3)
title("exactextract", line = -4)
invisible(NULL)
}
par(mar = rep(1.8, 4))
layout(matrix(c(1, 2), ncol = 2))
sp <- wk::as_wkt(silicate::inlandwaters)

fun(sp, 33, .1)

The parameter n is the number of pixels we want in the plane. So, first let's do a sensible default, a raster that is 10x10.

What do we see here? The yellow is good old terra::rasterize(), a single pixel raster in which was detected our two polygons

i <- 29
i <- i - 1
par(bg = "grey70")
fun(5)

![Image](https://github.com/user-attachments/assets/805bfbcd-b07b-42f9-8750-2fe6dc8b524b)