--- title: "bigcurve and s2: drawing the sphere on a flat map" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{bigcurve and s2: drawing the sphere on a flat map} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, 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)) ``` ```{r 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 ``` If you already use s2 (directly, or through sf, where it is the default for longitude and latitude data) you may wonder where bigcurve fits. The short answer: **they agree on what a line between two points on the Earth is, and they answer different questions about it.** * s2 answers *"what is this shape on the sphere, and how does it relate to other shapes?"* Does it intersect, what is its area, what is the union, how far apart are these? * bigcurve answers *"how many points do I need to draw this shape in a given map projection so that it looks right, and no more?"* In every figure below s2 is teal, bigcurve is black, and the planar reading of the same coordinates is a grey dashed line. One line runs through the whole guide: along latitude 60, from longitude 0 to 100. ## At a glance | | s2 | bigcurve | |---|---|---| | Job | geometry on the sphere: predicates, areas, unions, distances | drawing: just enough vertices for one map | | Edges | great circles | great circles | | Works in | 3D unit vectors, no projection | the plane of the target map | | Tolerance measured | on the sphere (metres converted to an angle) | on the map, in the map's units | | Projections for tessellating | plate carree, Mercator, orthographic | anything PROJ can do | | Output | s2 geography, or projected coordinates | longitude, latitude, so any tool can project it | | Input | s2 geography | matrices, segment meshes, wkpool topology (so any wk-readable data) | | Planar edges to geodesic | yes, on import (`planar = TRUE`) | no | ```{r setup} library(bigcurve) library(s2) the_line <- rbind(c(0, 60), c(100, 60)) line <- s2_make_line(the_line[, 1], the_line[, 2]) ``` ## What is a line between two points? Give a computer two points in longitude and latitude and ask for the line between them. There are two common answers. 1. The **planar** answer: draw a straight line in longitude, latitude space, as if the coordinates were x and y on graph paper. This was what sf and GEOS did for years, and what most file formats quietly assume. 2. The **geodesic** answer: follow the shortest path over the surface of the Earth, a great circle arc. s2 takes the second answer for every edge, and so does bigcurve. The planar middle of our line is (50, 60). The great circle middle is a long way north of it: ```{r 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) ``` ```{r 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) ``` More than a thousand kilometres. Neither answer is wrong: they are two different shapes that happen to share their vertices. What matters is knowing which one your data means. s2 makes the geodesic reading the default; bigcurve assumes it. ## Does bigcurve draw s2's line? Yes. bigcurve inserts its new vertices exactly on s2's arc. Densify the line and ask s2 how far each new vertex is from it: ```{r 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)) ``` A few nanometres, which is floating point noise. The added vertices change nothing about the shape on the sphere; they only make it explicit for a flat drawing. ## Where is the work done? s2 does all its work on the sphere itself. Internally every point is a unit vector in three dimensions, and there is no map projection anywhere in sight. That is why s2 has no trouble at the poles or at the antimeridian: on a sphere there is no edge of the map to fall off. bigcurve does its work in the **plane of the map you are going to draw**. For each edge it projects the great circle midpoint with PROJ and asks how far it lands from the straight chord between the projected endpoints. If it is further than `tolerance` (in the map's own units, usually metres) the edge is split and both halves are checked again. This is the adaptive resampling that Mike Bostock wrote for d3-geo. So the same edge gets **different** numbers of vertices in different projections, because the curve is a different shape on each map: ```{r 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))) ``` ## Where do the vertices go? s2 is not blind to this problem. Its C++ library has an *edge tessellator* that turns geodesic edges into chains of straight edges in a projection, and the R package exposes it when it writes data out, for plate carree (longitude, latitude), Mercator and orthographic maps. Here is our line tessellated for Mercator by both packages, each asked for 1 km: ```{r 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)) ``` ```{r 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) ``` In Mercator the two come out close, in count and in placement: both are solving the same problem in the same spirit. The difference is in what "1 km" means. ## What does a tolerance mean? s2 measures its tolerance **on the sphere**; bigcurve measures it **on the map**. On a map that stretches distances, as Mercator does by a factor of two at latitude 60, those are not the same thing. And s2 can only aim at its three projections: for any other map the best it can do is tessellate in longitude and latitude and hope that is fine enough. The figure measures, for every drawn chord, how far it strays on the map from the true projected line. ```{r 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) ``` On the left, s2's 1 km on the sphere becomes up to `r sprintf("%.1f", max(e_s2_merc$err) / 1000)` km on the Mercator map near the ends of the line, while bigcurve stays under 1 km everywhere. On the right, s2's longitude, latitude tessellation happens to be far finer than this map needs: `r nrow(v_s2_ll)` vertices for errors of tens of metres, where `r nrow(v_bc_st)` keep bigcurve under 1 km. On another map, or another line, the same tessellation could just as easily be too coarse. bigcurve's tolerance is always "how far may the drawn line stray on *this* map". The "planar edges to geodesic" row of the table deserves a sentence too. If your data really are planar in longitude and latitude (a coastline digitised on a plate carree map, say), s2 can densify on the way *in* (`s2_geog_from_wkb(..., planar = TRUE)`) so that its geodesic edges follow the planar ones. bigcurve never does this. It always assumes the edges you hand it are great circles. ## Using them together Because both packages agree on what an edge is, a natural pipeline is s2 for the analysis and bigcurve for the drawing. Here is a "box" written with four corners, from longitude 100 to 160 and latitude -35 to -20. We cut Australia with it in s2, then densify what the cut produced for a map in plain longitude and latitude: ```{r 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))) ``` ```{r 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) ``` Read as s2 reads it, the box's north side is not latitude -20: it is a great circle that dips to `r sprintf("%.1f", min(gc_path(c(100, -20), c(160, -20))[, 2]))` in the middle, and that is the line s2 cut along. The coastline is already finely detailed and is left alone; the new vertices go only along the long sides the box contributed. Drawing the grey reading while analysing the black one is how a map ends up disagreeing with its own analysis. To draw the result in any other projection, turn the pool back into geometry and project it. The fast general path is wk and PROJ: build the transform once with `PROJ::proj_trans_create()` and stream the geometry through `wk::wk_transform()`. Any other tool that projects (sf, terra) works too, since the densified vertices are ordinary longitude, latitude points. Remember to densify for the map you are drawing: ```{r 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) ``` ## What else is the same, and what is not? * **Both are spherical.** s2 uses a sphere of radius about 6371 km; so does bigcurve's midpoint. Neither follows the ellipsoid. For drawing a map the difference is small (the WGS84 geodesic from Hobart to London runs up to about 23 km to one side of the great circle). For distances and paths on the ellipsoid use the geographiclib package (GeographicLib is also inside sf and geosphere). * **Neither clips.** If part of a line is beyond the horizon of an orthographic map, s2 can intersect with `s2_hemisphere()` first; bigcurve leaves unprojectable edges alone and never errors. * **The antimeridian.** s2 has no seam at all. bigcurve keeps whatever longitudes you gave it: an edge along 180 stays at 180, an edge written from 170 to -170 is refined through 180 continuously. If you need the two copies of the seam to be one boundary, densify first and then fuse with `wkpool::merge_coincident(period_x = 360)`. * **Operations.** Union, intersection, buffer, distance, area and containment are s2's job. bigcurve has one verb, `densify()`, and it does not change the shape. ## In one sentence Use s2 to decide *what* the shape is; use bigcurve to draw *that* shape, honestly, in whatever projection your map needs.