An Interface to Rust's 'geo' Library.
rsgeo
rsgeo is an interface to the Rust libraries geo-types and geo. geo-types implements pure rust geometry primitives. The geo library adds additional algorithm functionalities on top of geo-types. This package lets you harness the speed, safety, and memory efficiency of these libraries. geo-types does not support Z or M dimensions. There is no support for CRS at this moment.
# install.packages(
# 'rsgeo',
# repos = c('https://josiahparry.r-universe.dev', 'https://cloud.r-project.org')
# )
library(rsgeo)
rsgeo works with vectors of geometries. When we compare this to sf this is always the geometry column which is a class sfc object (simple feature column).
# get geometry from sf
data(guerry, package = "sfdep")
polys <- guerry[["geometry"]] |>
sf::st_cast("POLYGON")
# cast to rust geo-types
rs_polys <- as_rsgeo(polys)
head(rs_polys)
#> <rs_POLYGON[6]>
#> [1] Polygon { exterior: LineString([Coord { x: 801150.0, y: 2092615.0 }, Coord...
#> [2] Polygon { exterior: LineString([Coord { x: 729326.0, y: 2521619.0 }, Coord...
#> [3] Polygon { exterior: LineString([Coord { x: 710830.0, y: 2137350.0 }, Coord...
#> [4] Polygon { exterior: LineString([Coord { x: 882701.0, y: 1920024.0 }, Coord...
#> [5] Polygon { exterior: LineString([Coord { x: 886504.0, y: 1922890.0 }, Coord...
#> [6] Polygon { exterior: LineString([Coord { x: 747008.0, y: 1925789.0 }, Coord...
Cast geometries to sf
sf::st_as_sfc(rs_polys)
#> Geometry set for 116 features
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 47680 ymin: 1703258 xmax: 1031401 ymax: 2677441
#> CRS: NA
#> First 5 geometries:
#> POLYGON ((801150 2092615, 800669 2093190, 80068...
#> POLYGON ((729326 2521619, 729320 2521230, 72928...
#> POLYGON ((710830 2137350, 711746 2136617, 71243...
#> POLYGON ((882701 1920024, 882408 1920733, 88177...
#> POLYGON ((886504 1922890, 885733 1922978, 88547...
Calculate the unsigned area of polygons.
bench::mark(
rust = unsigned_area(rs_polys),
sf = sf::st_area(polys),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 rust 55.6µs 57.65µs 16411. 3.8KB 0
#> 2 sf 1.36ms 1.44ms 649. 786.9KB 8.42
Find centroids
bench::mark(
centroids(rs_polys),
sf::st_centroid(polys),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 centroids(rs_polys) 174.95µs 213µs 3720. 3.8KB 9.53
#> 2 sf::st_centroid(polys) 2.43ms 2.6ms 359. 892.9KB 4.70
Extract points coordinates
coords(rs_polys) |>
head()
#> x y line_id polygon_id
#> 1 801150 2092615 1 1
#> 2 800669 2093190 1 1
#> 3 800688 2095430 1 1
#> 4 800780 2095795 1 1
#> 5 800589 2096112 1 1
#> 6 800333 2097190 1 1
Plot the polygons and their centroids
plot(rs_polys)
plot(centroids(rs_polys), add = TRUE)
Calculate a distance matrix. Note that there is often floating point error differences so check = FALSE in this case.
pnts <- centroids(rs_polys)
pnts_sf <- sf::st_as_sfc(pnts)
bench::mark(
rust = distance_euclidean_matrix(pnts, pnts),
sf = sf::st_distance(pnts_sf, pnts_sf),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 rust 323.53µs 573.06µs 1540. 108KB 4.08
#> 2 sf 3.48ms 3.69ms 256. 351KB 0
Simplify geometries.
x <- rs_polys
x_simple <- simplify_geoms(x, 5000)
plot(x_simple)
bench::mark(
rust = simplify_geoms(rs_polys, 500),
sf = sf::st_simplify(polys, FALSE, 500),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 rust 6.29ms 6.76ms 141. 1.91KB 0
#> 2 sf 8.52ms 9.02ms 108. 1.24MB 2.08
Union geometries with union_geoms(). Some things sf is better at! One of which is performing unary unions of complex geometries.
plot(union_geoms(rs_polys))
bench::mark(
union_geoms(rs_polys),
sf::st_union(polys),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 union_geoms(rs_polys) 205ms 209ms 4.78 0B 0
#> 2 sf::st_union(polys) 120ms 134ms 7.49 921KB 0
We can cast between geometries as well.
lns <- cast_geoms(rs_polys, "linestring")
Some unions are faster when using rsgeo vectors like linestrings.
lns_sf <- sf::st_cast(polys, "LINESTRING")
bench::mark(
union_geoms(lns),
sf::st_union(lns_sf),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 union_geoms(lns) 117.5µs 174µs 4275. 0B 0
#> 2 sf::st_union(lns_sf) 87.8ms 94ms 10.7 2.46MB 2.68
Find the closest point to a geometry
close_pnt <- closest_point(
rs_polys,
geom_point(800000, 2090000)
)
plot(rs_polys[1])
plot(close_pnt, pch = 15, add = TRUE)
Find the haversine destination of a point, bearing, and distance. Compare to the very fast geosphere destination point function.
bench::mark(
rust = haversine_destination(geom_point(10, 10), 45, 10000),
Cpp = geosphere::destPoint(c(10, 10), 45, 10000),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 rust 5.86µs 7.34µs 120442. 3.2KB 12.0
#> 2 Cpp 17.06µs 19.11µs 34150. 11.8MB 34.2
origin <- geom_point(10, 10)
destination <- haversine_destination(origin, 45, 10000)
plot(c(origin, destination), col = c("red", "blue"))
Find intermediate point on a great circle.
middle <- haversine_intermediate(origin, destination, 1/2)
plot(origin)
plot(destination, add = TRUE, col = "red")
plot(middle, add = TRUE, col = "blue")
Find extreme coordinates with extreme_coords()
france <- union_geoms(rs_polys)
plot(france)
plot(extreme_coords(france)[[1]], add = TRUE, pch = 15)
Get bounding rectangles
rects <- bounding_rect(rs_polys)
plot(rects)
Convex hulls
convex_hull(rs_polys) |>
plot()
Expand into constituent geometries as a list of geometry vectors
expand_geoms(rs_polys) |>
head()
#> [[1]]
#> <rs_LINESTRING[1]>
#> [1] LineString([Coord { x: 801150.0, y: 2092615.0 }, Coord { x: 800669.0, y: 2...
#>
#> [[2]]
#> <rs_LINESTRING[2]>
#> [1] LineString([Coord { x: 729326.0, y: 2521619.0 }, Coord { x: 729320.0, y: 2...
#> [2] LineString([Coord { x: 647667.0, y: 2468296.0 }, Coord { x: 647777.0, y: 2...
#>
#> [[3]]
#> <rs_LINESTRING[1]>
#> [1] LineString([Coord { x: 710830.0, y: 2137350.0 }, Coord { x: 711746.0, y: 2...
#>
#> [[4]]
#> <rs_LINESTRING[1]>
#> [1] LineString([Coord { x: 882701.0, y: 1920024.0 }, Coord { x: 882408.0, y: 1...
#>
#> [[5]]
#> <rs_LINESTRING[1]>
#> [1] LineString([Coord { x: 886504.0, y: 1922890.0 }, Coord { x: 885733.0, y: 1...
#>
#> [[6]]
#> <rs_LINESTRING[1]>
#> [1] LineString([Coord { x: 747008.0, y: 1925789.0 }, Coord { x: 746630.0, y: 1...
We can flatten the resultant geometries into a single vector using flatten_geoms()
expand_geoms(rs_polys) |>
flatten_geoms() |>
head()
#> <rs_LINESTRING[6]>
#> [1] LineString([Coord { x: 801150.0, y: 2092615.0 }, Coord { x: 800669.0, y: 2...
#> [2] LineString([Coord { x: 729326.0, y: 2521619.0 }, Coord { x: 729320.0, y: 2...
#> [3] LineString([Coord { x: 647667.0, y: 2468296.0 }, Coord { x: 647777.0, y: 2...
#> [4] LineString([Coord { x: 710830.0, y: 2137350.0 }, Coord { x: 711746.0, y: 2...
#> [5] LineString([Coord { x: 882701.0, y: 1920024.0 }, Coord { x: 882408.0, y: 1...
#> [6] LineString([Coord { x: 886504.0, y: 1922890.0 }, Coord { x: 885733.0, y: 1...
Combine geometries into a single multi- geometry
combine_geoms(lns)
#> <rs_LINESTRING[1]>
#> [1] MultiLineString([LineString([Coord { x: 801150.0, y: 2092615.0 }, Coord { ...
Spatial predicates
x <- rs_polys[1:5]
intersects_sparse(x, rs_polys)
#> [[1]]
#> [1] 1 48 50 92 94
#>
#> [[2]]
#> [1] 2 7 63 78 80 81 98 101
#>
#> [[3]]
#> [1] 3 20 27 53 77 84 94
#>
#> [[4]]
#> [1] 4 5 30 107 109
#>
#> [[5]]
#> [1] 4 5 30 48
Notes
Right now plotting is done using wk by first casting the rsgeo into an sfc object.