## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 6
)

## ----setup--------------------------------------------------------------------
library(sf)
library(terra)
library(weightedVoronoi)

## -----------------------------------------------------------------------------
# Terrain-anisotropy example (self-contained)
crs_use2 <- 32636

boundary_sf2 <- st_sf(
  id = 1,
  geometry = st_sfc(
    st_polygon(list(rbind(
      c(0, 0),
      c(1200, 0),
      c(1200, 900),
      c(800, 900),
      c(800, 500),
      c(400, 500),
      c(400, 900),
      c(0, 900),
      c(0, 0)
    ))),
    crs = crs_use2
  )
)

dem_aniso <- terra::rast(
  ext = terra::ext(terra::vect(boundary_sf2)),
  resolution = 20,
  crs = terra::crs(terra::vect(boundary_sf2))
)

xy2 <- terra::crds(dem_aniso, df = TRUE)
terra::values(dem_aniso) <- xy2$x * 10

pts2 <- st_sf(
  village = c("A", "B"),
  population = c(1, 1),
  geometry = st_sfc(
    st_point(c(200, 450)),
    st_point(c(950, 450))
  ),
  crs = crs_use2
)

out_geo_iso2 <- weighted_voronoi_domain(
  points_sf = pts2,
  weight_col = "population",
  boundary_sf = boundary_sf2,
  res = 20,
  distance = "geodesic",
  dem_rast = dem_aniso,
  use_tobler = TRUE,
  anisotropy = "none",
  verbose = FALSE
)

out_geo_aniso2 <- weighted_voronoi_domain(
  points_sf = pts2,
  weight_col = "population",
  boundary_sf = boundary_sf2,
  res = 20,
  distance = "geodesic",
  dem_rast = dem_aniso,
  use_tobler = TRUE,
  anisotropy = "terrain",
  uphill_factor = 3,
  downhill_factor = 1.2,
  verbose = FALSE
)

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))

plot(st_geometry(boundary_sf2), border = "black", lwd = 2, main = "Isotropic DEM geodesic")
plot(out_geo_iso2$polygons["generator_id"], add = TRUE)
plot(st_geometry(pts2), add = TRUE, pch = 21, bg = "red")

plot(st_geometry(boundary_sf2), border = "black", lwd = 2, main = "Terrain-anisotropic geodesic")
plot(out_geo_aniso2$polygons["generator_id"], add = TRUE)
plot(st_geometry(pts2), add = TRUE, pch = 21, bg = "red")

par(oldpar)

## -----------------------------------------------------------------------------
crs_u <- "EPSG:32636"

boundary_u <- st_sf(
  geometry = st_sfc(st_polygon(list(rbind(
    c(0, 0),
    c(200, 0),
    c(200, 200),
    c(0, 200),
    c(0, 0)
  )))),
  crs = crs_u
)

points_u <- st_sf(
  population = c(0.01, 0.02),
  geometry = st_sfc(
    st_point(c(60, 100)),
    st_point(c(140, 100))
  ),
  crs = crs_u
)

out_u <- weighted_voronoi_uncertainty(
  points_sf = points_u,
  weight_col = "population",
  boundary_sf = boundary_u,
  n_sim = 30,
  weight_sd = 0.8,
  distance = "geodesic",
  geodesic_engine = "multisource",
  weight_model = "additive",
  verbose = FALSE,
  seed = 1
)

plot(out_u$entropy, main = "Entropy")
plot(terra::vect(points_u), add = TRUE, pch = 21, col = "black", bg = "red")

## -----------------------------------------------------------------------------
pts_t1 <- points_u
pts_t2 <- points_u
pts_t2$population <- c(0.02, 0.01)

out_time <- weighted_voronoi_time(
  points_list = list(time1 = pts_t1, time2 = pts_t2),
  weight_col = "population",
  boundary_sf = boundary_u,
  distance = "geodesic",
  geodesic_engine = "multisource",
  weight_model = "additive",
  verbose = FALSE
)

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))
plot(out_time$change_map_first_last, main = "Change map")
plot(out_time$persistence, main = "Persistence")
par(oldpar)