usage example

[mdsumner]

library(textures)
quad2segment <- function(mesh) {
  ib <- mesh$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 = "-")))
  rgl::mesh3d(vertices = mesh$vb, segments = allseg[, !dup])
}

prop_globe <- function(x) {
  ## 360 * 180
  prod(diff(x)[c(1, 3)])/64800
}

EXT <- c(-180, 180, -65, -30)
#segmesh <- quad2segment(quad(c(360, 180)/15, extent = c(-180, 180, -75, 75)))
segmesh <- quad2segment(quad(c(18, 12), extent = EXT))

l <- vector("list", ncol(segmesh$is))
prj <- "+proj=laea +lat_0=0 +lon_0=147"
for (i in seq_along(l)) {
  e <- try(bigcurve:::bisect(t(segmesh$vb[1:2, segmesh$is[,i] ]), prj, dist = sqrt(prop_globe(EXT) * pi * 6378137 )/8))
  if (inherits(e, "try-error")) next; 

  l[[i]] <- rbind(e, NA)
}

x <- terra::project(do.call(rbind, l), to  = prj, from = "OGC:CRS84")
dim(x)
plot(x, pch = ".", asp = 1)
lines(x)
points(x, cex = 0.2, pch = 19)
m <- do.call(cbind, maps::map(xlim = EXT[1:2], ylim = EXT[3:4], plot = FALSE)[1:2])
lines(terra::project(m, to = prj, from  = "OGC:CRS84"))

motivation is the graticule package, and the various pieces are linked in this issue: https://github.com/hypertidy/textures/issues/12

and Dewey's example (a good one for describing why we want loxodromes and when)

#POLYGON ((-70 50, 0 50, 0 70, -35 52, -70 70, -70 50))

EXT <- c(-70, 0, 50, 70)
segmesh <- rgl::mesh3d(cbind(c(-70, 0, 0, -35, -70), 
                  c(50, 50, 70, 52, 70), 0, 1), segments = c(1, 2, 2, 3, 3, 4, 4, 5, 5, 1))

l <- vector("list", ncol(segmesh$is))
prj <- "+proj=ortho +lon_0=-35 +lat_0=58"
for (i in seq_along(l)) {
  e <- try(bigcurve:::bisect(t(segmesh$vb[1:2, segmesh$is[,i] ]), prj, dist = sqrt(prop_globe(EXT) * pi * 6378137 )/16))
  if (inherits(e, "try-error")) next; 

  l[[i]] <- rbind(e, NA)
}

x <- terra::project(do.call(rbind, l), to  = prj, from = "OGC:CRS84")
dim(x)
plot(x, pch = ".", asp = 1)
lines(x)
points(x, cex = 0.2, pch = 19)
m <- do.call(cbind, maps::map(xlim = EXT[1:2], ylim = EXT[3:4], plot = FALSE)[1:2])
lines(terra::project(m, to = prj, from  = "OGC:CRS84"))

[mdsumner]

things have moved on a lot, there's an internal "grat()" for something that might replace {graticule}, this example just compares the points added by s2 (assuming a tesselation value so you get uniform density) with those added by bigcurve, because the lines needed more vertices to represet them in a given projection.

I haven't unpacked s2 yet so it's hard to hack the projections.

  library(s2)
lons = seq(-0, 140, by = 10); 
lats = seq(-60, -40, by = 5)
proj <-"+proj=laea +lat_0=-50"
plot(grct <- wk::as_rct(wk::grd(wk::rct(min(lons), min(lats), max(lons), max(lats)), dx = 10, dy = 5)))

s2_ortho <- function(x) {
  centre <- as.numeric(s2:::last_plot_env$centre); 
  sprintf("+proj=ortho +lon_0=%f +lat_0=%f +R=6378137", centre[1], centre[2])
}

s2_proj <- function(x) {
  terra::project(as.matrix(x), to = s2_ortho(), from = "OGC:CRS84")/6378137
}

s2_plot(
 geog <-  s2_geog_from_wkb(
    wk::as_wkb(grct),
    planar = TRUE,
    tessellate_tol_m = 1e3
  ),
  border = "red"
)

s2_plot(s2_data_countries(), add = TRUE)
points(s2_proj(wk::wk_coords(geog)[, c("x", "y")]))
g <- bigcurve:::grat(lons = seq(-0, 140, by = 10), lats = seq(-60, -40, by = 10), proj = "+proj=laea +lat_0=-50 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

#> [1] 204956.8 749263.0
#> [1] 291010424691 400398718039

plot(g)
points(wk::wk_coords(g)[, c("x", "y")])

^{Created on 2023-06-19 with reprex v2.0.2}