Skip to content

Instantly share code, notes, and snippets.

@mdsumner

mdsumner/voropoly.md

Created August 28, 2024 11:14
Show Gist options
  • Save mdsumner/8db5ac01e47fa86f10e7ebc372e0ebda to your computer and use it in GitHub Desktop.
Save mdsumner/8db5ac01e47fa86f10e7ebc372e0ebda to your computer and use it in GitHub Desktop.
Download ZIP
Raw
voropoly.md 
sfx <- silicate::inlandwaters[6, ]
## devtools::install_github("hypertidy/anglr") ## in desperate need of reboot, alignment with wk/geos etc. 
ps <- anglr::as_pslg(sfx)
tr <- triangulate(ps)
plot(tr$VP, pch = ".")
plot(tr, add = TRUE)
points(tr$VP, pch = "+")
E <- tr$VE[tr$VE[,2] > 0, ]
lines(tr$VP[t(cbind(E, NA)), ], col = "red", lwd = 2)
@mdsumner
Copy link
Author

mdsumner commented Aug 29, 2024 •
edited
Loading 

    library(geos)
    library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.9.0, PROJ 9.3.1; sf_use_s2() is TRUE
    library(RTriangle)
    library(anglr)
    library(terra)
#> terra 1.7.81
    sfx <- silicate::minimal_mesh[1, ]
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
## devtools::install_github("hypertidy/anglr") ## in desperate need of reboot, alignment with wk/geos etc. 
ps <- anglr::as_pslg(sfx)
tr <- RTriangle::triangulate(ps)

plot(vpoly <-  geos::geos_voronoi_polygons(sfx), lwd = 5)
plot(sfx, col = scales::alpha("red", 0.2), add = TRUE)

points(tr$VP, pch = "+")
E <- tr$VE[tr$VE[,2] > 0, ]
lines(tr$VP[t(cbind(E, NA)), ], col = "red", lwd = 2)

## and now use terra intersect which does what you expect

tv <- terra::vect(sf::st_as_sf(vpoly))
ts <- terra::vect(sfx)
plot(tv)

plot(ts)

plot(intersect(tv, ts))

Created on 2024-08-29 with reprex v2.0.2

@mdsumner
Copy link
Author

all with terra

    ## do it all with terra
    library(terra)
#> terra 1.7.81
    f <- system.file("ex/lux.shp", package="terra")
    v <- vect(f)
    vo <- voronoi(v)
    
        
plot(intersect(v, vo))

plot(intersect(aggregate(v), voronoi(aggregate(v))))

Created on 2024-08-29 with reprex v2.0.2

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment