R/addContour.R

In cysouw/qlcVisualize: Visualization for Quantitative Language Comparison

addContour <- function(heights, points, window, crs, add = TRUE, polygons = FALSE,
                       levels = c(.4, .45, .5), grid = 5e4, ...) {

  # prepare data
  data <- sf::st_sf(geometry = sf::st_geometry(points), feature = heights, crs = crs)
  w <- sf::as_Spatial(sf::st_geometry(window))

  # make grid inside polygons
  grd <- as.data.frame(sp::spsample(w, "regular", n = grid))
  names(grd) <- c("X", "Y")
  grd <- sf::st_as_sf(grd, coords = c("X", "Y"), crs = crs)
  grd <- stars::st_rasterize(grd)

  # single polygon gets zero outside instead of NA
  grd$ID[grd$ID == 0] <- NA

  # ordinary kriging
  m <- automap::autofitVariogram(feature ~ 1, data)
  k <- gstat::krige(feature ~ 1, data, newdata = grd, model = m$var_model)

  if (add) {
    # draw contour with default lwd
    z <- k["var1.pred",,]
    if (is.null(list(...)$lwd)) {
      lwd <- rev(seq(2, 0.5, length.out = length(levels)))
      contour(z, add = TRUE, drawlabels = FALSE,
              levels = levels, lwd = lwd, ...)
    } else {
      contour(z, add = TRUE, drawlabels = FALSE,
              levels = levels, ...)
    }
  } else {
    # return linestrings or polygons
    cont <- stars::st_contour(k["var1.pred"], na.rm = FALSE,
                            contour_lines = !polygons, breaks = levels)
    names(cont)[1] <- "levels"
    return(cont)
  }
}