explore hausdorf distances

hausdorff-explorer.R

library(sf)

#' Plots vertices where Hausdorff distance is calculated from
#' 
#' @param x any single geometry
#' @param y any single geometry
#' @returns two points as an `sfc` object that are the locations where Hausdorff distance
#' is calculated from
plot_hd_points <- function(x, y) {
  
  # if they're not sfc convert them 
  if (!inherits(x, "sfc")) x <- st_as_sfc(x)
  if (!inherits(x, "sfc")) y <- st_as_sfc(y)
  
  x_pnts <- st_cast(x, "POINT")
  y_pnts <- st_cast(y, "POINT")
  
  d_mat <- st_distance(x_pnts, y_pnts)
  
  d <- as.numeric(st_distance(x, y, which = "Hausdorff"))
  
  origin_pnt_index <- which(d_mat == d, arr.ind = TRUE)
  
  pnts <- c(x_pnts[origin_pnt_index[1]], y_pnts[origin_pnt_index[2]])
  
  plot(c(x, y))
  plot(pnts, add = TRUE, col = "red", pch = 16, cex = 2)
}

# use some data from sf
# you can use your own shapeful path here if you'd like
nc <- st_read(system.file("shape/nc.shp", package = "sf")) |> 
  st_geometry() |> 
  # remove crs to not mess with unit conversions
  st_set_crs(NA)

plot_hd_points(nc[1], nc[2])
plot_hd_points(nc[1], nc[18])


# use data from sfdep package
geo <- sfdep::guerry$geometry

plot_hd_points(geo[1], geo[37])
plot_hd_points(geo[1], geo[36])
plot_hd_points(geo[1], geo[67])
plot_hd_points(geo[1], geo[69])