Package {potentiomap}

Title: Build Potentiometric Surfaces and Flow Arrows
Version: 0.1.0
Author: Elvin Cordero [aut, cre]
Maintainer: Elvin Cordero <elvin.cordero@seamountgeo.com>
Description: Builds potentiometric surface products from groundwater monitoring data. The package prepares groundwater observations from direct water-level measurements or depth-to-water data paired with land-surface elevations, interpolates thin-plate spline surfaces by default, supports alternative and user-supplied interpolation methods, exports raster and contour products, and derives hydraulic-gradient flow arrows. Raster operations use methods from Hijmans (2025) <doi:10.32614/CRAN.package.terra>, thin-plate spline interpolation uses methods from Nychka et al. (2021) <doi:10.5065/D6W957CT>, and geostatistical interpolation uses methods from Pebesma (2004) <doi:10.1016/j.cageo.2004.03.012>.
License: GPL-3
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.3.3
Depends: R (≥ 4.1)
Imports: fields, grDevices, graphics, gstat, sf, stats, terra
Suggests: testthat (≥ 3.0.0)
Config/testthat/edition: 3
NeedsCompilation: no
Packaged: 2026-05-25 04:22:47 UTC; ec
Repository: CRAN
Date/Publication: 2026-05-29 10:10:02 UTC

Extract arrow base or tip points

Description

Extract arrow base or tip points

Usage

ps_arrow_vertices(arrows, which = c("last", "first"), out_file = NULL)

Arguments

arrows

A line SpatVector or path to a line vector file.

which

"first" for arrow bases or "last" for arrow tips.

out_file

Optional output vector path.

Value

A point SpatVector.

Examples

data("synthetic_wells")
pts <- ps_make_points(synthetic_wells, "x", "y", "gw_elevation",
                      "well_id", "EPSG:26916")
s <- ps_interpolate(pts, methods = "IDW", grid_res = 100)
arrows <- ps_flow_arrows(s$IDW, res_factor = 4, scale = 60)
tips <- ps_arrow_vertices(arrows$arrows, which = "last")
tips

Create contours from a surface raster

Description

Create contours from a surface raster

Usage

ps_contours(surface, interval = 1, levels = NULL)

Arguments

surface

A terra::SpatRaster potentiometric surface.

interval

Contour interval in map elevation units.

levels

Optional explicit contour levels. When supplied, interval is ignored.

Value

A line terra::SpatVector of contours.

Examples

data("synthetic_wells")
pts <- ps_make_points(synthetic_wells, "x", "y", "gw_elevation",
                      "well_id", "EPSG:26916")
s <- ps_interpolate(pts, methods = "IDW", grid_res = 100)
ctr <- ps_contours(s$IDW, interval = 1)
ctr

Export surfaces, contours, and quicklook PNGs

Description

Export surfaces, contours, and quicklook PNGs

Usage

ps_export_surfaces(
  surfaces,
  out_dir,
  out_stub = "gw",
  contour_interval = 1,
  points = NULL,
  write_raster = TRUE,
  write_contours = TRUE,
  write_png = TRUE
)

Arguments

surfaces

A named list of SpatRaster objects, such as the result of ps_interpolate().

out_dir

Output directory.

out_stub

File prefix.

contour_interval

Contour interval.

points

Optional observation points to draw on quicklook figures.

write_raster, write_contours, write_png

Choose which outputs to write.

Value

A data frame listing written files.

Examples

data("synthetic_wells")
pts <- ps_make_points(synthetic_wells, "x", "y", "gw_elevation",
                      "well_id", "EPSG:26916")
s <- ps_interpolate(pts, methods = "IDW", grid_res = 100)
out <- ps_export_surfaces(s, points = pts, out_dir = tempdir())
out

Generate hydraulic-gradient flow arrows

Description

Derives slope, aspect, hydraulic gradient, and downgradient arrows from a potentiometric surface raster.

Usage

ps_flow_arrows(
  surface,
  res_factor = 7,
  scale = 50,
  min_gradient = 1e-05,
  log_gradient = FALSE,
  log_arrow = FALSE,
  out_dir = NULL,
  out_stub = "gw"
)

Arguments

surface

A groundwater elevation SpatRaster.

res_factor

Factor used to thin arrows by resampling to a coarser grid.

scale

Arrow length multiplier.

min_gradient

Gradients below this value are dropped.

log_gradient

Store log1p() transformed gradient in the output raster.

log_arrow

Use log1p() transformed gradient for arrow lengths.

out_dir

Optional output directory. When supplied, files are written.

out_stub

File prefix used when writing outputs.

Value

A list with raster, points, and arrows.

Examples

data("synthetic_wells")
pts <- ps_make_points(synthetic_wells, "x", "y", "gw_elevation",
                      "well_id", "EPSG:26916")
s <- ps_interpolate(pts, methods = "IDW", grid_res = 100)
arrows <- ps_flow_arrows(s$IDW, res_factor = 4, scale = 60)
arrows$arrows

Interpolate potentiometric surfaces

Description

Creates one raster per requested interpolation method. The default method is thin-plate spline ("TPS"). Other built-in methods are inverse distance weighting ("IDW"), ordinary kriging ("OK"), and universal kriging with quadratic drift ("UK"). Advanced users can also pass named custom interpolation functions through custom_methods.

Usage

ps_interpolate(
  points,
  value = "Z",
  methods = "TPS",
  grid_res = NULL,
  template = NULL,
  mask = NULL,
  padding = NULL,
  idw_power = 2,
  idw_nmax = 15,
  tps_lambda = NULL,
  kr_auto_cutoff = TRUE,
  kr_cutoff = NA_real_,
  kr_width = NA_real_,
  custom_methods = NULL,
  x = "x",
  y = "y",
  name_col = NULL,
  crs = NULL
)

Arguments

points

A point SpatVector, sf object, or coordinate table with a data column to interpolate.

value

Data column name when points is not already standardized. Defaults to "Z".

methods

Character vector of interpolation methods. Built-in values are "TPS", "IDW", "OK", and "UK". Names supplied in custom_methods can also be used.

grid_res

Output raster cell size in map units.

template

Optional template SpatRaster; overrides grid_res, padding, and mask extent construction.

mask

Optional AOI polygon used to crop and mask output rasters.

padding

Padding added around the convex hull extent when building a template from points.

idw_power, idw_nmax

IDW power and maximum neighbors.

tps_lambda

Thin-plate spline smoothing parameter. NULL lets fields::Tps() choose by GCV.

kr_auto_cutoff

Use automatic variogram cutoff and lag width.

kr_cutoff, kr_width

Manual variogram cutoff and lag width.

custom_methods

Optional named list of custom interpolation functions. Each function is called as fun(points, template, grid) and must return either a SpatRaster matching template or a numeric vector with one value per template cell.

x, y, name_col, crs

Used when points is a coordinate table.

Value

A named list of SpatRaster surfaces.

Examples

data("synthetic_wells")
pts <- ps_make_points(synthetic_wells, "x", "y", "gw_elevation",
                      "well_id", "EPSG:26916")
surfaces <- ps_interpolate(pts, grid_res = 100)
names(surfaces)

Make groundwater observation points

Description

Convert a coordinate table, sf point object, or terra vector to a SpatVector with standard Z and Name fields.

Usage

ps_make_points(data, x = "x", y = "y", value, name_col = NULL, crs = NULL)

Arguments

data

A data frame, sf object, or terra::SpatVector.

x, y

Coordinate column names for tabular data.

value

Groundwater elevation column name.

name_col

Optional well or station name column.

crs

Coordinate reference system for tabular data, such as "EPSG:26916".

Value

A point terra::SpatVector with standardized attributes.

Examples

data("synthetic_wells")
pts <- ps_make_points(
  synthetic_wells,
  x = "x", y = "y",
  value = "gw_elevation",
  name_col = "well_id",
  crs = "EPSG:26916"
)
pts

Build potentiometric points from depth-to-water measurements

Description

Calculates groundwater elevation as surface elevation minus depth to water. Surface elevation can come from a DEM raster, a column in the depth table, or separate surface-elevation points. Separate surface points are matched by name when possible; otherwise their elevations are interpolated to the depth points with inverse distance weighting.

Usage

ps_potentiometric_points(
  data,
  x = "x",
  y = "y",
  depth_col,
  surface = NULL,
  surface_col = NULL,
  name_col = NULL,
  surface_name_col = name_col,
  crs = NULL,
  idw_power = 2
)

Arguments

data

Depth-to-water observations as a data frame, sf, or terra::SpatVector.

x, y

Coordinate column names for tabular depth data.

depth_col

Depth-to-water column name. Positive values are assumed to be depth below land surface.

surface

A DEM SpatRaster, separate surface-elevation observations, or NULL when surface_col is used.

surface_col

Surface-elevation column in data, or in surface when surface is a point/table object.

name_col

Optional name column in data.

surface_name_col

Optional name column in separate surface observations.

crs

CRS for tabular depth data.

idw_power

Power used when interpolating separate surface points.

Value

A point terra::SpatVector with surface_elevation, depth_to_water, and Z groundwater elevation fields.

Examples

data("synthetic_wells")
data("synthetic_dem")
gw <- ps_potentiometric_points(
  synthetic_wells,
  x = "x", y = "y",
  depth_col = "depth_to_water",
  surface = synthetic_dem,
  name_col = "well_id",
  crs = "EPSG:26916"
)
head(terra::values(gw))

Draw a quicklook surface plot

Description

Draw a quicklook surface plot

Usage

ps_quicklook(
  surface,
  contours = NULL,
  points = NULL,
  file = NULL,
  title = "Potentiometric surface",
  label_points = TRUE,
  width = 1600,
  height = 1200,
  res = 180
)

Arguments

surface

A SpatRaster.

contours

Optional contour SpatVector.

points

Optional observation point SpatVector.

file

Optional PNG output path. When NULL, plots to the active device.

title

Plot title.

label_points

Label points with Name and Z.

width, height, res

PNG dimensions and resolution.

Value

Invisibly returns file.

Examples

data("synthetic_wells")
pts <- ps_make_points(synthetic_wells, "x", "y", "gw_elevation",
                      "well_id", "EPSG:26916")
s <- ps_interpolate(pts, methods = "IDW", grid_res = 100)
ps_quicklook(s$IDW, points = pts, title = "Synthetic IDW")

Make the sample area-of-interest polygon

Description

Make the sample area-of-interest polygon

Usage

ps_sample_aoi()

Value

A terra::SpatVector polygon in EPSG:26916.

Examples

aoi <- ps_sample_aoi()
aoi

Smooth a potentiometric surface raster

Description

Applies a focal moving-window smoother to a potentiometric surface raster. This can be useful when an interpolated surface is technically valid but too locally rough for contour development or hydraulic-gradient visualization.

Usage

ps_smooth_surface(
  surface,
  window_size = 3,
  method = c("mean", "median"),
  weights = NULL,
  iterations = 1,
  na.rm = TRUE,
  preserve_na = TRUE,
  filename = "",
  overwrite = FALSE
)

Arguments

surface

A terra::SpatRaster potentiometric surface.

window_size

Odd integer window size used when weights is NULL.

method

Smoothing statistic. Supported values are "mean" and "median".

weights

Optional odd-dimension numeric matrix of focal weights. NA values in the matrix are ignored by terra::focal().

iterations

Number of smoothing passes.

na.rm

Ignore missing values inside the focal window.

preserve_na

Preserve the original NA footprint after smoothing.

filename

Optional output raster filename.

overwrite

Overwrite filename when it exists.

Value

A smoothed terra::SpatRaster.

Examples

data("synthetic_wells")
pts <- ps_make_points(synthetic_wells, "x", "y", "gw_elevation",
                      "well_id", "EPSG:26916")
s <- ps_interpolate(pts, grid_res = 100)
smoothed <- ps_smooth_surface(s$TPS, window_size = 5)
smoothed

Synthetic DEM raster

Description

A small artificial DEM matching synthetic_wells, stored as a packed terra raster. Use terra::rast(synthetic_dem) to unpack it.

Usage

synthetic_dem

Format

A terra::PackedSpatRaster with one layer named surface_elevation.

Examples

data("synthetic_dem")
dem <- terra::rast(synthetic_dem)
dem

Synthetic surface elevation measurement points

Description

Artificial land-surface elevation points for demonstrating workflows that do not start with a DEM raster.

Usage

synthetic_surface_points

Format

A data frame with coordinate, surface-elevation, and name columns.

Examples

data("synthetic_surface_points")
head(synthetic_surface_points)

Synthetic groundwater monitoring wells

Description

A real-looking artificial monitoring dataset with coordinates, land-surface elevation, depth to water, and calculated groundwater elevation.

Usage

synthetic_wells

Format

A data frame with 32 rows and 6 columns:

well_id

Synthetic well identifier.

x

Easting in EPSG:26916 map units.

y

Northing in EPSG:26916 map units.

surface_elevation

Land-surface elevation.

depth_to_water

Depth to groundwater below land surface.

gw_elevation

Groundwater elevation.

Examples

data("synthetic_wells")
head(synthetic_wells)