rebuild with wk and bigcurve

[mdsumner]

more examples here: https://gist.github.com/mdsumner/e6487e3272d11d13860e00ce24033361

## helper to just project verts
mkproj <- function(crs, source = "OGC:CRS84") {
  
  function(x) {
    reproj::reproj_xy(x, crs, source = source)
  }
}
## convert quad index (efficient mesh index/vert form of raster) to unique segments
quad2segment <- function(ib) {
  
  allseg <- cbind(ib[1:2, ], ib[2:3, ], ib[3:4, ], ib[c(4, 1), ])
  dup <- duplicated(pmax(paste(allseg[1, ], allseg[2, ], sep = "-"), 
                         paste(allseg[2, ], allseg[1, ], sep = "-")))
  allseg[, !dup]
}

#pak::pak("hypertidy/reproj")
#pak::pak("hypertidy/quad")
#pak::pak("hypertidy/bigcurve")
#pak::pak("maps")
#pak::pak("scales")
library(quad)
library(bigcurve)
dm <- c(18, 9)
ib <- matrix(quad_index(dm, ydown = TRUE), nrow = 4L) + 1
vb <- matrix(quad_vert(dm, ydown = TRUE), nrow = 2L)
## scale to world
xlim <- c(50, 180)
ylim <-  c(-75, 15)
vb[1, ] <- scales::rescale(vb[1, ], xlim)
vb[2, ] <- scales::rescale(vb[2, ],ylim)

## get map data
m <- do.call(cbind, maps::map(xlim = xlim, ylim = ylim, plot = F)[1:2])
plot(t(vb), pch = ".")
polygon(t(vb[, rbind(ib, ib[1, ], NA)]))


is <- quad2segment(ib)

## we should cull segments that span the entire domain
crs <- "+proj=omerc +lonc=140 +lat_0=-42 +alpha=45"
proj <- mkproj(crs)

## cull the segments for projected length
vp <- proj(t(vb))
## a bad segment is longer than earth radius-ish (we could quant this to get a heuristic)
dst <- apply(is, 2, function(.x) dist(vp[.x, ]))
bad <- (dst > 6378137 * 3.5) 
is <- is[,!bad]  
wk0 <- vector("list", ncol(is))
#plot(proj(t(vb)), asp = 1)
for (i in 1:ncol(is)) {
  bs <- bigcurve:::bisect(t(vb[,is[, i]]), crs, 5e3)
  wk0[[i]] <- wk::wk_polygon(wk::xy(bs[,1], bs[,2]))
  #points(proj(bs), pch = 19, cex = .2)
  #lines(proj(bs), pch = 19, cex = .2)
}

plot(wk0 <- do.call(c, wk0))
points(wk::wk_coords(wk0)[, c("x", "y"), drop = FALSE])