Skip to content

Instantly share code, notes, and snippets.

@JosiahParry
Created June 18, 2022 01:38
Show Gist options
  • Save JosiahParry/03509cdc6e7dbbc5590c1884c8e215cd to your computer and use it in GitHub Desktop.
Save JosiahParry/03509cdc6e7dbbc5590c1884c8e215cd to your computer and use it in GitHub Desktop.
development of alpha hull algorithm
nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"))
geo <- sf::st_centroid(nc) |>
sf::st_geometry()
tris <- geo |>
geos::geos_make_collection() |>
geos::geos_delaunay_triangles()
# this is what i would've preferred (and now do)
x <- geos::geos_unnest(tris)
# Identify the circumscribed circle
circumscribe_circle <- function(geom) {
tri_xy <- sf::st_coordinates(sf::st_as_sf(geom))[1:3,1:2]
circum_circ <- tripack::circumcircle(tri_xy[,1], tri_xy[,2])
circum_circ
}
# apply the above function to all triangles
all_circs <- do.call("rbind", lapply(x, circumscribe_circle))
# create new wk circles object
circs <- wk::new_wk_crc(
x = list(x = unlist(all_circs$x),
y = unlist(all_circs$y),
r = unlist(all_circs$r)),
crs = attr(x, "crs")
)
# pull the radii from the circles
radii <- unclass(circs)[["r"]]
# sort
radii_sorted <- sort(radii, TRUE)
# calculate alpha
alpha_og <- 1 / radii_sorted[1]
# plot tesselations
plot(x[radii < alpha_og])
@rsbivand
Copy link

Right, alphahull is a good choice. However, I suspect it assumes planar coordinates:

geo <- sf::st_centroid(nc) |> sf::st_transform("EPSG:32119") |>
  sfheaders::sf_to_df() |> dplyr::select(x, y)
a <- alphahull::ahull(geo[,1], geo[,2], 1000000)
plot(a, asp=1)

where changing the tension on alpha gives varying effects. Please also note that tripack should be replaced by interp or deldir as it has an encumbered license.

@JosiahParry
Copy link
Author

@rsbivand & @mdsumner thank you! I've found that using sf::st_jitter()enables alphahull to work. This is after either using unique() on a matrix, or distinct on a data frame.

@mdsumner
Copy link

Roger please ... if course it does, it's not rocket science

use a projection, we're just getting basic shit to work the actual reality is immaterial at this level, where it's pure geometry

@mdsumner
Copy link

of course jitter would work, you'll get unique vertices which is exactly what the message says is required

don't get sucked into thinking the longlat thing has anything to do with it, or any of the other specious sf bs

@JosiahParry
Copy link
Author

Let's not get worked up over triangles and circles. :) I've tried both long lat and projected coords, both work after jittering.

@mdsumner
Copy link

I'm sorry Jos, you know it's just numbers ... I'm not going to participate in any nonsense about it being "non-planar" (when that, while real, is irrelevant)

I don't have any patience for that idiocy and that's why I'll just dip out 🙏

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