## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 4.5, dpi = 96,
                      out.width = "100%")
have <- all(vapply(c("s2", "wk", "wkpool", "maps"), requireNamespace,
                   logical(1), quietly = TRUE))
knitr::opts_chunk$set(eval = have)
## the last drawing uses the PROJ package; only run it where PROJ can
## find its database and build a transform
proj_ok <- have && requireNamespace("PROJ", quietly = TRUE) &&
  isTRUE(tryCatch(!is.null(PROJ::proj_trans_create("OGC:CRS84", "EPSG:3857")),
                  error = function(e) FALSE))

## ----style, include = FALSE---------------------------------------------------
## one colour per package in every figure
col_s2 <- "#1a8a8a"
col_bc <- "black"
col_planar <- "grey50"
col_coast <- "grey78"
R <- 6378137
coast <- do.call(cbind, maps::map(plot = FALSE)[1:2])
coast[which(coast[, 1] > 180), ] <- NA

## ----setup--------------------------------------------------------------------
library(bigcurve)
library(s2)
the_line <- rbind(c(0, 60), c(100, 60))
line <- s2_make_line(the_line[, 1], the_line[, 2])

## ----midpoint-----------------------------------------------------------------
s2_as_text(s2_interpolate_normalized(line, 0.5))

## how far is the planar middle from the true line, in metres?
s2_distance(s2_lnglat(50, 60), line)

## ----fig-meanings, echo = FALSE, fig.width = 9, fig.height = 4.6--------------
gc_path <- function(a, b, n = 400) {
  l <- s2_make_line(c(a[1], b[1]), c(a[2], b[2]))
  p <- s2_interpolate_normalized(l, seq(0, 1, length.out = n))
  as.matrix(wk::wk_coords(p)[, c("x", "y")])
}
true_ll <- gc_path(the_line[1, ], the_line[2, ])
planar_ll <- cbind(seq(0, 100, length.out = 200), 60)
mid_planar <- c(50, 60)
mid_gc <- as.numeric(wk::wk_coords(s2_interpolate_normalized(line, 0.5))[, c("x", "y")])
gap_km <- s2_distance(s2_lnglat(mid_planar[1], mid_planar[2]),
                      s2_lnglat(mid_gc[1], mid_gc[2])) / 1000

op <- par(mfrow = c(1, 2), mar = c(0.5, 0.5, 2.2, 0.5))
plot(NA, xlim = c(-20, 120), ylim = c(35, 85), asp = 1 / cos(60 * pi / 180),
     axes = FALSE, xlab = "", ylab = "")
lines(coast, col = col_coast)
lines(planar_ll, col = col_planar, lty = 2, lwd = 2)
lines(true_ll, col = col_bc, lwd = 2)
segments(mid_planar[1], mid_planar[2], mid_gc[1], mid_gc[2], col = "firebrick", lwd = 1.5)
points(rbind(mid_planar, mid_gc), pch = 21, bg = "white", cex = 1.1)
points(the_line, pch = 19)
text(mean(c(mid_gc[1], mid_planar[1])), mean(c(mid_gc[2], mid_planar[2])),
     sprintf("%.0f km", gap_km), pos = 4, col = "firebrick")
title("Longitude and latitude as x and y", font.main = 1, cex.main = 1)

ortho <- "+proj=ortho +lon_0=50 +lat_0=55 +R=6378137"
ang <- seq(0, 2 * pi, length.out = 361)
plot(R * cos(ang), R * sin(ang), type = "n", asp = 1, axes = FALSE, xlab = "", ylab = "")
polygon(R * cos(ang), R * sin(ang), col = "grey97", border = "grey60")
lines(rproj_xy(coast, ortho), col = col_coast)
lines(rproj_xy(planar_ll, ortho), col = col_planar, lty = 2, lwd = 2)
lines(rproj_xy(true_ll, ortho), col = col_bc, lwd = 2)
mids <- rproj_xy(rbind(mid_planar, mid_gc), ortho)
segments(mids[1, 1], mids[1, 2], mids[2, 1], mids[2, 2], col = "firebrick", lwd = 1.5)
points(mids, pch = 21, bg = "white", cex = 1.1)
points(rproj_xy(the_line, ortho), pch = 19)
title("The same two points on the globe", font.main = 1, cex.main = 1)
legend("bottomleft", bty = "n", cex = 0.85,
       lty = c(2, 1), lwd = 2, col = c(col_planar, col_bc),
       legend = c("planar: straight in lon, lat", "geodesic: great circle (s2, bigcurve)"))
par(op)

## ----on-the-arc---------------------------------------------------------------
d <- densify(the_line, "+proj=laea +lon_0=50 +lat_0=60", tolerance = 1000)
nrow(d)
max(s2_distance(s2_lnglat(d[, 1], d[, 2]), line))

## ----per-projection-----------------------------------------------------------
c(laea  = nrow(densify(the_line, "+proj=laea +lon_0=50 +lat_0=60", tolerance = 1000)),
  stere = nrow(densify(the_line, "+proj=stere +lat_0=90", tolerance = 1000)),
  merc  = nrow(densify(the_line, "+proj=merc", tolerance = 1000)))

## ----tessellate---------------------------------------------------------------
merc <- "+proj=merc +R=6378137"
tol <- 1000

## s2, in its own Mercator, then back to longitude and latitude
m <- wk::wk_handle(line, wk::wkb_writer(),
                   s2_projection = s2_projection_mercator(),
                   s2_tessellate_tol = tol / s2_earth_radius_meters())
v_s2 <- rproj_xy(as.matrix(wk::wk_coords(wk::wkb(m))[, c("x", "y")]),
                 "OGC:CRS84", source = merc)

## bigcurve, for the same Mercator
v_bc <- densify(the_line, merc, tolerance = tol)

c(s2 = nrow(v_s2), bigcurve = nrow(v_bc))

## ----fig-where, echo = FALSE, fig.width = 9, fig.height = 5.2-----------------
along <- function(ll) {
  as.numeric(s2_distance(s2_lnglat(the_line[1, 1], the_line[1, 2]),
                         s2_lnglat(ll[, 1], ll[, 2]))) / s2_length(line)
}
layout(matrix(1:2, 2), heights = c(3, 1.4))
op <- par(mar = c(0.5, 0.5, 2.2, 0.5))
xy_true <- rproj_xy(true_ll, merc)
pad <- 6e5
plot(NA, xlim = range(xy_true[, 1]) + c(-pad, pad),
     ylim = range(xy_true[, 2]) + c(-pad, pad), asp = 1,
     axes = FALSE, xlab = "", ylab = "")
lines(rproj_xy(coast, merc), col = col_coast)
lines(xy_true, col = "grey60", lwd = 4)
points(sweep(rproj_xy(v_s2, merc), 2, c(0, 1.2e5), "+"), pch = 19, cex = 0.45, col = col_s2)
points(sweep(rproj_xy(v_bc, merc), 2, c(0, 1.2e5), "-"), pch = 19, cex = 0.45, col = col_bc)
title(sprintf("One line in Mercator, 1 km tolerance: s2 %d vertices, bigcurve %d",
              nrow(v_s2), nrow(v_bc)), font.main = 1, cex.main = 1)
legend("bottomright", bty = "n", pch = 19, col = c(col_s2, col_bc), cex = 0.85,
       legend = c("s2 tessellation (drawn just above the line)",
                  "bigcurve (drawn just below)"))
par(mar = c(3.2, 7.5, 0.5, 0.5))
plot(NA, xlim = c(0, 1), ylim = c(0.5, 2.5), axes = FALSE, xlab = "", ylab = "")
segments(along(v_s2), 1.65, along(v_s2), 2.35, col = col_s2)
segments(along(v_bc), 0.65, along(v_bc), 1.35, col = col_bc)
axis(1, at = seq(0, 1, 0.25), labels = c("start", "", "half way", "", "end"))
axis(2, at = c(2, 1), labels = c("s2", "bigcurve"), las = 1, tick = FALSE)
mtext("position of each vertex along the line", side = 1, line = 2.1, cex = 0.85)
par(op)

## ----fig-tolerance, echo = FALSE, fig.width = 9, fig.height = 4.2-------------
map_error <- function(ll, prj) {
  t <- along(ll)
  xy <- rproj_xy(ll, prj)
  out <- numeric(nrow(ll) - 1L)
  for (i in seq_along(out)) {
    arc <- rproj_xy(gc_path(ll[i, ], ll[i + 1, ], n = 64), prj)
    a <- xy[i, ]; b <- xy[i + 1, ]
    d <- b - a; L <- sqrt(sum(d^2))
    p <- sweep(arc, 2, a)
    out[i] <- max(abs(d[1] * p[, 2] - d[2] * p[, 1]) / L)
  }
  data.frame(t0 = t[-length(t)], t1 = t[-1], err = out)
}
step_lines <- function(e, ...) {
  segments(e$t0, pmax(e$err, 1), e$t1, pmax(e$err, 1), ...)
}
stere <- "+proj=stere +lat_0=90 +lon_0=50 +R=6378137"
v_s2_ll <- as.matrix(wk::wk_coords(wk::wkt(
  s2_as_text(line, planar = TRUE, tessellate_tol_m = tol)))[, c("x", "y")])
v_bc_st <- densify(the_line, stere, tolerance = tol)
e_s2_merc <- map_error(v_s2, merc)
e_bc_merc <- map_error(v_bc, merc)
e_s2_st <- map_error(v_s2_ll, stere)
e_bc_st <- map_error(v_bc_st, stere)

op <- par(mfrow = c(1, 2), mar = c(3.4, 4.2, 3.2, 0.8))
ylim <- c(1, max(c(e_s2_merc$err, e_bc_merc$err, e_s2_st$err, e_bc_st$err)) * 1.5)
panel <- function(es2, ebc, main, s2_label) {
  plot(NA, xlim = c(0, 1), ylim = ylim, log = "y", axes = FALSE,
       xlab = "", ylab = "")
  abline(h = tol, col = "grey60", lty = 3)
  step_lines(ebc, col = col_bc, lwd = 3.5)
  step_lines(es2, col = col_s2, lwd = 1.8)
  axis(1, at = seq(0, 1, 0.5), labels = c("start", "half way", "end"))
  at <- c(1, 10, 100, 1000, 1e4, 1e5)
  at <- at[at <= ylim[2]]
  axis(2, at = at, labels = c("1 m", "10 m", "100 m", "1 km", "10 km", "100 km")[seq_along(at)], las = 1)
  mtext("drawn chord to true line, on the map", side = 2, line = 3.2, cex = 0.8)
  title(main, font.main = 1, cex.main = 1)
  legend("bottomleft", bty = "n", lwd = 2, col = c(col_s2, col_bc), cex = 0.8,
         legend = c(sprintf("s2: %s (%d vertices)", s2_label, nrow(es2) + 1L),
                    sprintf("bigcurve: 1 km on the map (%d vertices)", nrow(ebc) + 1L)))
}
panel(e_s2_merc, e_bc_merc, "Mercator (s2 can tessellate here)",
      "1 km on the sphere")
panel(e_s2_st, e_bc_st, "North polar stereographic (s2 cannot)",
      "1 km in lon, lat, then projected")
par(op)

## ----together-----------------------------------------------------------------
box <- s2_make_polygon(c(100, 160, 160, 100), c(-35, -35, -20, -20))
aus <- s2_data_countries("Australia")
cut <- s2_intersection(aus, box)

pool <- wkpool::establish_topology(wk::as_wkb(cut))
out <- densify(pool, "+proj=eqc")
c(before = nrow(wkpool::pool_vertices(pool)),
  after  = nrow(wkpool::pool_vertices(out)))

## ----fig-cut, echo = FALSE, fig.width = 9, fig.height = 6---------------------
v0 <- wkpool::pool_vertices(pool)
v1 <- wkpool::pool_vertices(out)
added <- v1[!v1$.vx %in% v0$.vx, ]
flat <- function(x) wk::wk_set_geodesic(wk::wk_set_crs(x, NULL), FALSE)
geom <- wkpool::cycles_to_wkb(wkpool::merge_coincident(out))

op <- par(mar = c(2.5, 3.2, 2.2, 0.5))
plot(NA, xlim = c(112, 155), ylim = c(-44, -10), asp = 1, xlab = "", ylab = "", las = 1)
plot(flat(wk::as_wkb(aus)), col = "grey95", border = col_coast, add = TRUE)
plot(flat(geom), col = "grey85", border = "grey40", add = TRUE)
abline(h = c(-20, -35), col = col_planar, lty = 2, lwd = 1.5)
text(c(154.5, 154.5), c(-20, -35), c("latitude -20", "latitude -35"),
     pos = 3, col = col_planar, cex = 0.85)
lines(gc_path(c(100, -20), c(160, -20)), col = col_bc, lwd = 1)
points(added$x, added$y, pch = 19, cex = 0.8, col = col_bc)
title(sprintf("Cutting Australia with the box: its north side is a great circle (%d vertices added)",
              nrow(added)), font.main = 1, cex.main = 0.95)
legend("bottomleft", bty = "n", cex = 0.85, lty = c(2, 1), pch = c(NA, 19),
       col = c(col_planar, col_bc), lwd = c(1.5, 1),
       legend = c("the box read as planar (parallels)",
                  "the box as s2 and bigcurve read it, with added vertices"))
par(op)

## ----together-draw, eval = proj_ok, fig.height = 3.5--------------------------
prj <- "+proj=laea +lon_0=130 +lat_0=-30"
geom <- wkpool::cycles_to_wkb(wkpool::merge_coincident(densify(pool, prj)))
drawn <- wk::wk_transform(geom, PROJ::proj_trans_create(wk::wk_crs(geom), prj))
dc <- wk::wk_coords(drawn)
op <- par(mar = rep(0.5, 4))
plot(range(dc$x), range(dc$y), type = "n", asp = 1,
     axes = FALSE, xlab = "", ylab = "")
plot(drawn, col = "grey85", add = TRUE)
par(op)