## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", fig.bg = "white",
  dev.args = list(bg = "white"), fig.width = 8, fig.height = 6,
  dpi = 144, fig.align = "center"
)

## ----prepare------------------------------------------------------------------
library(potentiomap)
data("synthetic_wells")

points <- ps_make_points(
  synthetic_wells,
  x = "x", y = "y", value = "gw_elevation", name_col = "well_id",
  crs = "EPSG:26916", head_unit = "m", output_unit = "m",
  vertical_datum = "synthetic example datum",
  surface_reference = "land_surface", metadata_mode = "strict"
)

## ----interpolate--------------------------------------------------------------
result <- ps_interpolate(
  points, methods = "IDW", grid_res = 300,
  return = "result", support = TRUE, support_max_distance = 1000
)
result
ps_diagnostics(result, "IDW")
result$support$summary

## ----first-map, fig.alt="Synthetic IDW potentiometric surface shaded blue at lower modeled heads and red at higher heads, with one-unit contours and observation wells."----
surface <- result$surfaces$IDW
contours <- ps_contours(surface, interval = 1)

par(bg = "white")
terra::plot(
  surface, col = hcl.colors(64, "RdYlBu", rev = TRUE),
  main = "First potentiometric surface"
)
terra::plot(contours, add = TRUE, col = "#263845", lwd = 1)
terra::plot(points, add = TRUE, pch = 21, bg = "white", cex = 0.75)

## ----export-------------------------------------------------------------------
output_dir <- file.path(tempdir(), "potentiomap-getting-started")
exports <- ps_export_surfaces(
  result, output_dir, out_stub = "synthetic", points = points,
  vector_format = "gpkg"
)
knitr::kable(data.frame(
  method = exports$method,
  raster = basename(exports$raster)
))