could be simpler, use the logic to remove the expand_xy step, just use vectors x, y and build indexes for those
## rect as stand in for "polygon" raster grid
## faster than expand.grid, for x,y case
expand_xy <- function (x, y) {
cbind(x, rep(y, each = length(x)))
}
## build every corner coordinate, the edges of each pixel
edges_xy <- function (dm, ...) {
xx <- seq(0, dm[1L], length = dm[1L] + 1L)
yy <- seq(dm[2L], 0, length = dm[2L] + 1L)
expand_xy(x = xx, y = yy)
}
## build the component indices (the xy grid is unique coordinates so we
## use grid logic to get interleaved corners xmin,xmax,ymin,ymax)
build_rct_index <- function(dm) {
nc <- dm[1L]
nr <- dm[2L]
xmin <- seq_len(nc) + rep(seq(0, length.out = nr, by = nc + 1), each = nc)
xmax <- xmin + 1
ymax <- xmax
ymin <- xmin + nc + 1
cbind(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax)
}
#sst <- raster::aggregate(raadtools::readsst(), fact = 5)
library(raster)
r <- raster::aggregate(anglr::gebco, fact = 16)
dm <- dim(r)[2:1]
#dm <- c(3, 4)
xy <- edges_xy(dm)
idx <- build_rct_index(dm)
## now say our grid has xlim, ylim (we started with c(0, ncol), c(0, nrow))
xlim <- c(-180, 180)
ylim <- c(-90, 90)
## because it's just a scaling, provides an opportunity for delaying
## this step, but we do it now for example
cxy <- cbind(scales::rescale(xy[,1], xlim),
scales::rescale(xy[,2], ylim))
library(wk)
g <- wk::rct(xmin = cxy[idx[,"xmin"], 1L], xmax = cxy[idx[,"xmax"], 1L],
ymin = cxy[idx[,"ymin"], 2L], ymax = cxy[idx[,"ymax"], 2L])
plot(g, col = palr::d_pal(raster::values(r)), border = NA)
maps::map(add = T)
could be simpler, use the logic to remove the expand_xy step, just use vectors x, y and build indexes for those