gebco <- "/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2024.tif"
library(gdalraster)
library(wk)
library(geos)
proj <- "+proj=spilhaus"
w <- warp(gebco, tf <- "/vsimem/spil.vrt", proj, cl_arg = c("-tr", 25000, 25000))
g <- wk::as_wkb( grd(bbox = rct(-180, -85, 180, 85), dx = 15, dy = 10, type = "polygons"), crs = "EPSG:4326")
## project the endpoints, create as linestrings and drop any that are too long
coords <- wk::wk_coords(geos::geos_densify(g, 2.5))
to_line_segs <- function(x) {
x <- as.matrix(x[, c("x", "y")])
i <- seq(1, nrow(x) - 1)
j <- seq(2, nrow(x))
txy <- rbind(c(t(cbind(x[i,1], x[j, 1]))),
c(t(cbind(x[i,2], x[j, 2]))))
wk::wk_linestring(wk::xy(txy[1, ], txy[2, ]), feature_id = rep(1:(nrow(x)-1), each = 2))
}
cxy <- gdalraster::transform_xy( coords[, c("x", "y")], srs_to = proj, srs_from = "EPSG:4326")
linesegslist <- lapply(split(tibble::tibble(x = cxy[,1], y = cxy[,2]), coords[, "feature_id"]), to_line_segs)
linesegs <- do.call(c, linesegslist)
bad <- geos::geos_length(linesegs) >= quantile(geos::geos_length(linesegs), .99) |
!(geos::geos_length(linesegs) > 0)
plot_raster(read_ds(new(GDALRaster, tf)))
## having some plotting problem I don't understand with wkb ...
plot(as_wkb(as_wkt(linesegs[!bad])), add = TRUE, col = "firebrick")
now with countries exploded to lines