cell2cell

[mdsumner]

cell2cell <- function(dm, blocksize, cell) {
  nc <- dm[1L]
  col0 <- (cell - 1L) %% nc          # 0-indexed col in source
  row0 <- (cell - 1L) %/% nc          # 0-indexed row in source
  
  tile_col <- col0 %/% blocksize[1L]  # 0-indexed tile col
  tile_row <- row0 %/% blocksize[2L]  # 0-indexed tile row
  
  n_tile_cols <- ceiling(dm[1L] / blocksize[1L])
  tile_row * n_tile_cols + tile_col + 1L
}

tile_lines <- function(dm, ex, blocksize, crs_src, crs_dst, n_dense = 100, reproj = TRUE) {
  nx <- dm[1L]; ny <- dm[2L]
  
  # tile boundary indices (0-indexed cell edges)
  x_breaks <- seq(ex[1], ex[2], length.out = nx + 1L)
  y_breaks <- seq(ex[3], ex[4], length.out = ny + 1L)
  
  tile_x <- x_breaks[seq(1L, nx + 1L, by = blocksize[1L])]
  tile_y <- y_breaks[seq(1L, ny + 1L, by = blocksize[2L])]
  
  # vertical lines (constant x), densified in y
  y_dense <- seq(ex[3], ex[4], length.out = n_dense)
  vlines <- lapply(tile_x, \(x) cbind(x = rep(x, n_dense), y = y_dense))
  
  # horizontal lines (constant y), densified in x
  x_dense <- seq(ex[1], ex[2], length.out = n_dense)
  hlines <- lapply(tile_y, \(y) cbind(x = x_dense, y = rep(y, n_dense)))
  
  all_lines <- c(vlines, hlines)
  if (!reproj) return(all_lines)
  # reproject each line
  invisible(lapply(all_lines, \(ln) reproj::reproj(ln, target = crs_dst, source = crs_src)))
}


src <- list(ex = c(-180, 180, -90, 90), dm = c(360, 180), crs  = "EPSG:4326")

dst <- list(ex = vaster::buffer_extent(reproj::reproj_extent(src$ex, tgt <- "+proj=laea", source = src$crs ), 1000), 
            dm = c(360, 360), crs = tgt)

library(terra)
r <- project(rast(sds::gebco()), rast(ext(src$ex), ncols = src$dm[1], nrows = src$dm[2], crs = src$crs), by_util = TRUE)
par(mfrow = c(2, 2))
plot(project(setValues(r, cell2cell(dim(r)[2:1], c(72, 36), seq_len(prod(dim(r)[2:1])))), 
             rast(ext(dst$ex), ncols = dst$dm[1], nrows = dst$dm[2], crs = dst$crs), by_util = TRUE))

lapply(tile_lines(src$dm, src$ex, c(72, 36), src$crs, dst$crs), lines)
plot(project(setValues(r, cell2cell(dim(r)[2:1], c(72, 36), seq_len(prod(dim(r)[2:1])))), 
             rast(ext(dst$ex), ncols = dst$dm[1], nrows = dst$dm[2], crs = dst$crs), by_util = TRUE))

lapply(tile_lines(dst$dm, dst$ex, c(72, 72), dst$crs, src$crs, reproj = FALSE), lines)

r2 <- project(rast(sds::gebco()), rast(ext(dst$ex), ncols = dst$dm[1], nrows = dst$dm[2], crs = dst$crs), by_util = TRUE)

plot(project(setValues(r2, cell2cell(dim(r2)[2:1], c(72, 72), seq_len(prod(dim(r2)[2:1])))), 
             rast(ext(src$ex), ncols = src$dm[1], nrows = src$dm[2], crs = src$crs), by_util = TRUE))

max_jump <- function(xy) {
  dx <- diff(xy[, 1]); dy <- diff(xy[, 2])
  max(sqrt(dx^2 + dy^2))
}

# cull lines where any single segment exceeds threshold
threshold <- diff(dst$ex[1:2]) * 0.3  # e.g. 30% of display width
reprojected_lines <- tile_lines(dst$dm, dst$ex, c(72, 72), dst$crs, src$crs)
good <- !(sapply(reprojected_lines, max_jump) < threshold)
good[is.na(good)] <- F
lapply(reprojected_lines[good], lines)
plot(project(setValues(r2, cell2cell(dim(r2)[2:1], c(72, 72), seq_len(prod(dim(r2)[2:1])))), 
             rast(ext(src$ex), ncols = src$dm[1], nrows = src$dm[2], crs = src$crs), by_util = TRUE))

lapply(tile_lines(src$dm, src$ex, c(72, 36), src$crs, dst$crs, reproj = FALSE), lines)