Package {bigcurve}

bigcurve: Adaptive Densification of Lines for Map Projections

Description

Densify lines adaptively, as required by the curvature a map projection imposes on them. Edges are read as great circle arcs, and vertices are inserted along them only where the projected curve deviates from a straight chord by more than a tolerance measured on the map, following the adaptive resampling method of the 'd3-geo' 'JavaScript' library described in Bostock (2016) https://github.com/d3/d3-geo. Projections are computed with the 'PROJ' coordinate transformation library (https://proj.org/), so any coordinate reference system it supports can be a target. Inputs can be coordinate matrices, segment meshes or vertex pool topology from 'wkpool', and edges shared by neighbouring polygons receive identical new vertices.

Author(s)

Maintainer: Michael D. Sumner mdsumner@gmail.com [copyright holder]

Authors:

• Michael D. Sumner mdsumner@gmail.com [copyright holder]

Other contributors:

• Simon Wotherspoon (developed the spherical track midpoint method that the great circle midpoint used here is based on) [contributor]

• Mike Bostock (author of the 'd3-geo' adaptive resampling method ported here; see inst/COPYRIGHTS) [contributor, copyright holder]

See Also

Useful links:

• https://github.com/hypertidy/bigcurve

• https://hypertidy.github.io/bigcurve/

• Report bugs at https://github.com/hypertidy/bigcurve/issues

Adaptively densify lines for a projection

Description

Insert vertices along great circle arcs exactly where a target projection bends them, so that straight chords drawn between output vertices stay within tolerance of the true projected curve. This is the adaptive resampling scheme of d3-geo (Bostock), generalized to paths and to segment meshes.

Usage

densify(
  x,
  target,
  source = NULL,
  tolerance = NULL,
  pixels = 2048L,
  max_depth = 16L,
  ellipsoid = FALSE
)

Arguments

Details

The error metric is evaluated in the plane of the target projection: the great circle midpoint of each candidate segment is forward-projected and compared to the segment's planar chord (perpendicular distance, plus a parametric check that catches strong shear). Segments are split recursively until they pass or max_depth is reached. No inverse transforms are used and each point is projected exactly once.

Coordinates are returned in the coordinates of source (geographic by default): project the result yourself for plotting, e.g. with rproj_xy() for a matrix or, for wk geometry (such as wkpool::segments_to_wkb() of a refined pool), with wk::wk_transform() and a PROJ::proj_trans_create() transform.

Points that cannot be projected (for example beyond the horizon of an orthographic projection) do not error: segments touching them are left un-densified, and refinement stops at unprojectable midpoints. Clipping is out of scope here.

Inserted vertices are emitted in the longitude frame of their input edge: an edge along longitude 180 refines at 180, its mirror at -180 refines at -180 (the two seams of a global layer remain distinct), and an edge written as crossing the antimeridian yields a locally continuous chain, possibly outside ⁠[-180, 180)⁠. Input vertices are always returned with their written values, so a crossing edge's path carries one numeric 360-degree step at its far endpoint (the same point on the sphere; PROJ is indifferent). Longitudes are never canonicalized; normalize afterwards if a canonical range is needed. Edges that cross the target projection's own discontinuity (which depends on lon_0) cannot form a valid planar chain and should be cut beforehand; refinement of such edges is bounded, not exploded, because arcs spanning less ground than tolerance are never split.

Every edge is a great circle arc on a sphere, the same reading of an edge that 's2' uses. Data whose edges are meant as straight lines in longitude and latitude describe a different shape, and densifying them here changes it.

The source must be geographic (longitude, latitude): reading projected coordinates (eastings, northings in metres) as great circle arcs is meaningless, so input whose coordinates fall outside lon ⁠[-360, 360]⁠ or lat ⁠[-90, 90]⁠ is rejected with an error. Reproject to lon,lat first. (Treating a projected source's edges as straight lines in its own plane is a separate model, not yet implemented.)

Shared edges stay shared. An edge refines to the same vertices whichever way it is written, bit for bit, so two neighbouring polygons get identical new vertices along their common boundary (and wkpool::merge_coincident() fuses them). An edge that appears twice with the same vertices, as in a merged pool or a segment mesh, is refined once and its new vertices are reused.

Value

the same kind of object as x, with vertices added: a matrix for a path, a list of matrices for a list, a mesh with appended vertices and refined is for a mesh, a wkpool with appended vertices (freshly minted .vx ids), refined directed segments, and .feature provenance carried onto every refined segment

References

https://bost.ocks.org/mike/example/

Examples

s <- segment(c(140, -164), c(-89, 80))
d <- densify(s, "+proj=laea +lon_0=147 +lat_0=-42")
nrow(d)
plot(rproj_xy(d, "+proj=laea +lon_0=147 +lat_0=-42"), type = "l", asp = 1)

## on the WGS84 ellipsoid rather than the sphere
if (requireNamespace("geographiclib", quietly = TRUE)) {
  e <- densify(d[c(1, nrow(d)), ], "+proj=laea +lon_0=147 +lat_0=-42",
               ellipsoid = TRUE)
  lines(rproj_xy(e, "+proj=laea +lon_0=147 +lat_0=-42"), col = "firebrick")
}

Projection string

Description

Lambert azimuthal equal area projection string, in old-style PROJ format.

Usage

laea(lon_0 = 0, lat_0 = 0, x_0 = 0, y_0 = 0)

Arguments

Details

Note that easting/northing are poor names for false X and Y offsets, but that is what they are called in PROJ.

Value

character string, projection

Examples

laea()
laea(147, -42)

Transform coordinates

Description

Forward-transform a two-column matrix of coordinates with PROJ. Authority codes are axis-normalized, so input is always x,y (lon,lat) order and output is easting,northing - even for 'EPSG:4326'.

Usage

rproj_xy(x, target, ..., source = NULL)

Arguments

Details

Points outside the domain of target are returned as NA rather than raising an error.

Value

two-column matrix of transformed coordinates

Examples

rproj_xy(cbind(147, -42), "+proj=laea +lon_0=147 +lat_0=-42")

Segment

Description

A two-point path as a two-column matrix, by default a random one.

Usage

segment(x = runif(2L, -180, 180), y = runif(2L, -90, 90))

Arguments

Value

two-column matrix

Examples

segment()
segment(c(0, 147), c(0, -42))