Geometric functions provided by the simple features package
Ivan Lizarazo
10.12.2023
This notebook illustrates several geometric functions provided by the simple features for R package.
More information about this package is available here.
1. Setup
library(tmap)
library(tmaptools)
library(sf)
library(readr)
library(ggplot2)2. Data
Let’s start reading the toy dataset created using geojson.io functions:
list.files()## [1] "007.gpx" "009.gpx"
## [3] "009TRACKS.gpkg" "9TRACKS.gpkg"
## [5] "campus.geojson" "LINES.gpkg"
## [7] "NUEVETRACKS.gpkg" "poligon.geojson"
## [9] "poligonos.gpkg" "POLYS.gpkg"
## [11] "puntos.geojson" "simple_features.html"
## [13] "simple_features.nb.html" "simple_features.Rmd"campus <- st_read("campus.geojson")## Reading layer `campus' from data source
## `/Users/ials/Documents/unal/GB2023/campus/campus.geojson' using driver `GeoJSON'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.09209 ymin: 4.632038 xmax: -74.07948 ymax: 4.645007
## Geodetic CRS: WGS 84poly <- st_read("poligon.geojson")## Reading layer `poligon' from data source
## `/Users/ials/Documents/unal/GB2023/campus/poligon.geojson'
## using driver `GeoJSON'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.09282 ymin: 4.635755 xmax: -74.08279 ymax: 4.644851
## Geodetic CRS: WGS 84ptos <- st_read("puntos.geojson")## Reading layer `puntos' from data source
## `/Users/ials/Documents/unal/GB2023/campus/puntos.geojson' using driver `GeoJSON'
## Simple feature collection with 17 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -74.09355 ymin: 4.633727 xmax: -74.08041 ymax: 4.645315
## Geodetic CRS: WGS 84campus$var1 <- sample(1:2, nrow(campus), replace = TRUE)campuspoly$var2 <- sample(1:2, nrow(poly), replace = TRUE)polyptos$var3 <- sample(1:8, nrow(ptos), replace=TRUE)ptos3. Plots
tm_shape(campus) +
tm_fill(col = "green") +
tm_shape(poly, is.master = FALSE) +
tm_fill(col="cyan") +
tm_layout(legend.outside = TRUE)
# 4. Checking geometries
Do two existing geometries intersect?
sf::st_intersects(campus,poly)## Sparse geometry binary predicate list of length 1, where the predicate
## was `intersects'
## 1: 1What is the intersection between two polygons?
inter1 = sf::st_intersection(campus,poly)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesPrint the intersection:
inter1Plot the intersection:
tm_shape(campus) +
tm_fill(col = "gray") +
tm_shape(inter1, is.master = FALSE) +
tm_fill(col="red") +
tm_layout(legend.outside = TRUE)
5. Getting intersection between polys and points
First, plot the geometries:
tm_shape(campus) +
tm_polygons('#f0f0f0f0', border.alpha = 0.2) +
tm_shape(ptos)+
tm_dots(alpha = 0.7, col='green', size=0.5)
Then, get the intersection geometry:
inter2 = sf::st_intersection(campus,ptos)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesThen, print the intersection geometry:
inter2Now, plot the intersection geometry:
tm_shape(campus) +
tm_polygons('#f0f0f0f0', border.alpha = 0.2) +
tm_shape(inter2) +
tm_dots(alpha = 0.7, col='red', size=0.5, title = "puntos dentro del campus")
6. More geometric functions
What is the “boundary” of the campus?
(bounds1 = st_boundary(campus))What is the “convex hull” of the campus?
(hull1 = st_convex_hull(campus))What is the “buffer”, with 50 meter radius, of the campus?
(buffer1 = st_buffer(campus, 50))Let’s plot the campus and the 50 m. buffer:
tm_shape(campus) +
tm_fill(col = "green") +
tm_shape(buffer1, is.master = FALSE) +
tm_fill(col="cyan", alpha=0.5) +
tm_layout(legend.outside = TRUE)
Now, plot the campus and the convex hull:
tm_shape(campus) +
tm_fill(col = "green") +
tm_shape(hull1, is.master=FALSE) +
tm_fill(col="orange", alpha=0.6) +
tm_layout(legend.outside = TRUE)
Then, plot the campus and its boundary:
tm_shape(campus) +
tm_fill(col = "green") +
tm_shape(bounds1, is.master=FALSE) +
tm_lines(col="red", ldw=2, scale=3) +
tm_layout(legend.outside = TRUE)
7. Geometry creation
Before trying a triangulation, we need to reproject the points located inside the campus:
st_crs(inter2)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]Let’s transform the points to the *IGAC - Origen Unico Nacional” CRS:
(new_ptos = st_transform(inter2, crs=9377)) Now, the triangulation task:
edges <- new_ptos |>
st_union() |>
st_triangulate(bOnlyEdges = TRUE) |>
st_cast("LINESTRING")edges## Geometry set for 15 features
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 4879297 ymin: 2070242 xmax: 4880210 ymax: 2071078
## Projected CRS: MAGNA-SIRGAS / Origen-Nacional
## First 5 geometries:## LINESTRING (4880007 2070959, 4880210 2070773)## LINESTRING (4879623 2071078, 4880007 2070959)## LINESTRING (4879297 2070543, 4879623 2071078)## LINESTRING (4879297 2070543, 4879406 2070406)## LINESTRING (4879406 2070406, 4879870 2070242)The plot:
tm_shape(campus) +
tm_fill(col = "gray") +
tm_shape(edges, is.master=FALSE) +
tm_lines(col="red", scale=2) +
tm_layout(legend.outside = TRUE)
Now, let’s try a Voronoi task:
voros <- new_ptos |>
st_union() |>
st_voronoi(bOnlyEdges = TRUE) |>
st_cast("LINESTRING")The plot:
tm_shape(campus) +
tm_fill(col = "gray") +
tm_shape(voros, is.master=FALSE) +
tm_lines(col="red", scale=2) +
tm_layout(legend.outside = TRUE)
7. Assignment
Your task is to try additional functions in your notebook. Use the Spatial manipulation with sf: :CHEAT SHEET as reference.
8. References
Edzer Pebesma, 2018. Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal 10:1, 439-446.
9. Reproducibility
sessionInfo()## R version 4.3.0 (2023-04-21)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.4.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.4.3 readr_2.1.4 sf_1.0-14 tmaptools_3.1-1
## [5] tmap_3.3-3
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.4 xfun_0.39 bslib_0.4.2 raster_3.6-20
## [5] htmlwidgets_1.6.2 lattice_0.21-8 tzdb_0.4.0 vctrs_0.6.3
## [9] tools_4.3.0 crosstalk_1.2.0 generics_0.1.3 parallel_4.3.0
## [13] tibble_3.2.1 proxy_0.4-27 fansi_1.0.4 highr_0.10
## [17] pkgconfig_2.0.3 KernSmooth_2.23-20 RColorBrewer_1.1-3 leaflet_2.1.2
## [21] lifecycle_1.0.3 compiler_4.3.0 munsell_0.5.0 terra_1.7-46
## [25] codetools_0.2-19 leafsync_0.1.0 stars_0.6-1 htmltools_0.5.5
## [29] class_7.3-21 sass_0.4.6 yaml_2.3.7 pillar_1.9.0
## [33] jquerylib_0.1.4 classInt_0.4-10 cachem_1.0.8 lwgeom_0.2-13
## [37] wk_0.8.0 abind_1.4-5 tidyselect_1.2.0 digest_0.6.31
## [41] dplyr_1.1.2 fastmap_1.1.1 grid_4.3.0 colorspace_2.1-0
## [45] cli_3.6.1 magrittr_2.0.3 base64enc_0.1-3 dichromat_2.0-0.1
## [49] XML_3.99-0.14 utf8_1.2.3 leafem_0.2.0 e1071_1.7-13
## [53] withr_2.5.0 scales_1.2.1 sp_2.0-0 rmarkdown_2.22
## [57] png_0.1-8 hms_1.1.3 evaluate_0.21 knitr_1.43
## [61] viridisLite_0.4.2 s2_1.1.4 rlang_1.1.1 Rcpp_1.0.11
## [65] glue_1.6.2 DBI_1.1.3 rstudioapi_0.15.0 jsonlite_1.8.7
## [69] R6_2.5.1 units_0.8-4