| Title: | 'LiDAR' Terrain Analysis and Change Detection from Above |
|---|---|
| Description: | Terrain change detection, cut and fill volume estimation, terrain profiling, flood inundation analysis, slope and aspect computation, hillshade generation, contour extraction, reclamation monitoring, erosion analysis, and engineering export (LandXML, STL) from 'LiDAR' (Light Detection and Ranging) point clouds and digital elevation models ('DEMs'). Applications include surface mine reclamation monitoring, sediment pond capacity tracking, highwall safety classification, and erosion channel detection. Built on 'lidR' for point cloud I/O and 'terra' for raster operations. Includes access utilities for 'KyFromAbove' cloud-native elevation data on Amazon Web Services ('AWS') <https://kyfromabove.ky.gov/>. Methods for terrain change detection and volume estimation follow Li and others (2005) <doi:10.1016/j.geomorph.2004.10.007>. |
| Authors: | Chris Lyons [aut, cre] |
| Maintainer: | Chris Lyons <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.0.0 |
| Built: | 2026-06-03 07:41:21 UTC |
| Source: | https://github.com/chrislyonsky/abover |
Identifies flat bench areas between highwall faces in surface mining terrain. Benches are detected as relatively flat areas (low slope) bounded by steep terrain (high slope), using slope segmentation and area filtering.
bench_detection(dem, max_slope = 10, min_area = 100, highwall_slope = 40)bench_detection(dem, max_slope = 10, min_area = 100, highwall_slope = 40)
dem |
A terra::SpatRaster representing the terrain surface. |
max_slope |
Numeric. Maximum slope in degrees for a cell to be
considered part of a bench. Default |
min_area |
Numeric. Minimum contiguous area (in map units
squared) for a bench. Smaller patches are removed. Default |
highwall_slope |
Numeric. Minimum slope for adjacent highwall
terrain. Default |
An sf::sf data frame of bench polygons with columns:
bench_id: sequential bench identifier
area_m2: bench area in square metres
mean_elev: mean elevation of the bench surface
mean_slope: mean slope within the bench
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) benches <- bench_detection(dem, max_slope = 15, min_area = 0) if (nrow(benches) > 0) plot(sf::st_geometry(benches))dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) benches <- bench_detection(dem, max_slope = 15, min_area = 0) if (nrow(benches) > 0) plot(sf::st_geometry(benches))
Converts a polygon boundary to a closed linestring and samples elevation values around the perimeter. Useful for inspecting highwall edges, dam crests, or property boundaries.
boundary_terrain_profile(dem, boundary, spacing = NULL)boundary_terrain_profile(dem, boundary, spacing = NULL)
dem |
A terra::SpatRaster representing the terrain surface. |
boundary |
An sf::sf polygon whose boundary (exterior ring) will be profiled. |
spacing |
Numeric. Distance between sample points, in CRS units
of |
A data frame with columns:
distance: distance along the boundary perimeter
elevation: sampled elevation value
x, y: coordinates of each sample point
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) boundary <- sf::st_read( system.file("extdata/boundary.gpkg", package = "aboveR"), quiet = TRUE ) bprof <- boundary_terrain_profile(dem, boundary) plot(bprof$distance, bprof$elevation, type = "l", xlab = "Perimeter Distance", ylab = "Elevation")dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) boundary <- sf::st_read( system.file("extdata/boundary.gpkg", package = "aboveR"), quiet = TRUE ) bprof <- boundary_terrain_profile(dem, boundary) plot(bprof$distance, bprof$elevation, type = "l", xlab = "Perimeter Distance", ylab = "Elevation")
Extracts change statistics for each polygon in a zone layer, computing cut volume, fill volume, net change, and descriptive statistics per zone.
change_by_zone(change_raster, zones, id_field)change_by_zone(change_raster, zones, id_field)
change_raster |
A terra::SpatRaster as returned by |
zones |
An sf::sf data frame of polygons defining analysis zones. |
id_field |
Character. Column name in |
An sf::sf data frame with columns:
zone identifier
cut_volume: total volume of material removed (m^3, positive)
fill_volume: total volume of material added (m^3, positive)
net_volume: fill - cut (m^3)
area_m2: zone area
mean_change: mean elevation change
max_cut: deepest cut (most negative value)
max_fill: highest fill (most positive value)
before <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) after <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) chg <- terrain_change(before, after) zones <- sf::st_read( system.file("extdata/zones.gpkg", package = "aboveR"), quiet = TRUE ) change_by_zone(chg, zones, id_field = "zone_id")before <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) after <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) chg <- terrain_change(before, after) zones <- sf::st_read( system.file("extdata/zones.gpkg", package = "aboveR"), quiet = TRUE ) change_by_zone(chg, zones, id_field = "zone_id")
Returns a blue-white-red color ramp centered at zero, suitable for terrain change visualisation. Cut (negative) is blue, fill (positive) is red.
change_colors(n = 100L)change_colors(n = 100L)
n |
Integer. Number of colors. Default |
Character vector of hex color values.
cols <- change_colors(10) image(matrix(1:10), col = cols)cols <- change_colors(10) image(matrix(1:10), col = cols)
Identifies steep terrain faces (highwalls) typical of surface mining operations by computing slope from a DEM and classifying cells that exceed a slope threshold. Returns a binary raster and optionally vectorised polygons of highwall zones.
classify_highwall(dem, slope_threshold = 60, min_area = 0, as_polygons = FALSE)classify_highwall(dem, slope_threshold = 60, min_area = 0, as_polygons = FALSE)
dem |
A terra::SpatRaster representing the terrain surface. |
slope_threshold |
Numeric. Minimum slope in degrees to classify
as highwall. Default |
min_area |
Numeric. Minimum contiguous area (in map units squared)
for a highwall zone. Smaller patches are removed. Default |
as_polygons |
Logical. Return vectorised polygons instead of a
raster? Default |
If as_polygons = FALSE, a terra::SpatRaster with values
1 (highwall) and NA (non-highwall). If as_polygons = TRUE, an
sf::sf data frame of highwall polygons with an area_m2 column.
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) hw <- classify_highwall(dem, slope_threshold = 5) terra::plot(hw)dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) hw <- classify_highwall(dem, slope_threshold = 5) terra::plot(hw)
Extracts contour lines at a specified interval and returns them as
an sf::sf object with an elevation attribute.
contour_lines(dem, interval = 5, levels = NULL)contour_lines(dem, interval = 5, levels = NULL)
dem |
A terra::SpatRaster representing the terrain surface. |
interval |
Numeric. Elevation interval between contour lines.
Default |
levels |
Numeric vector. Specific elevations to contour. If
provided, |
An sf::sf data frame with LINESTRING geometry and an
elevation column.
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) cl <- contour_lines(dem, interval = 2) plot(sf::st_geometry(cl))dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) cl <- contour_lines(dem, interval = 2) plot(sf::st_geometry(cl))
Identifies potential erosion channels (rills, gullies) by computing flow accumulation and filtering cells where accumulated flow exceeds a threshold. Optionally returns vectorised channel lines.
detect_channels(dem, threshold = 100, as_lines = FALSE)detect_channels(dem, threshold = 100, as_lines = FALSE)
dem |
|
threshold |
Numeric. Minimum flow accumulation value to classify
as a channel. Default |
as_lines |
Logical. Return channel centrelines as sf::sf lines?
Default |
If as_lines = FALSE, a terra::SpatRaster with channel
cells = 1, other = NA. If as_lines = TRUE, an sf::sf LINESTRING
object of detected channels.
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) channels <- detect_channels(dem, threshold = 50) terra::plot(channels)dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) channels <- detect_channels(dem, threshold = 50) terra::plot(channels)
Computes the volume of material between a surface DEM and a
reference surface (e.g., a design grade or pre-mining DEM),
optionally clipped to a boundary polygon. Reports cut and fill
volumes separately and documents the integration method used.
estimate_volume( surface, reference, boundary = NULL, method = c("trapezoidal", "simpson") )estimate_volume( surface, reference, boundary = NULL, method = c("trapezoidal", "simpson") )
surface |
A terra::SpatRaster representing the actual surface. |
reference |
A terra::SpatRaster or single numeric value representing the reference elevation. If numeric, a constant reference plane at that elevation is used. |
boundary |
An optional sf::sf polygon to clip both surfaces to before computation. |
method |
Character. Volume integration method: |
A list with components:
cut_volume_m3: volume of material removed (positive)
fill_volume_m3: volume of material added (positive)
net_volume_m3: fill minus cut
area_m2: total analysed area
mean_depth_m: mean absolute difference
max_cut_m: deepest cut
max_fill_m: highest fill
method: integration method used
surface <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) reference <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) boundary <- sf::st_read( system.file("extdata/boundary.gpkg", package = "aboveR"), quiet = TRUE ) vol <- estimate_volume(surface, reference, boundary) cat("Net volume:", vol$net_volume_m3, "m3\n")surface <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) reference <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) boundary <- sf::st_read( system.file("extdata/boundary.gpkg", package = "aboveR"), quiet = TRUE ) vol <- estimate_volume(surface, reference, boundary) cat("Net volume:", vol$net_volume_m3, "m3\n")
Converts a DEM raster to a LandXML file containing a TIN (Triangulated Irregular Network) surface. LandXML is widely supported by civil engineering software including Autodesk Civil 3D, Bentley OpenRoads, and Trimble Business Center.
export_landxml( dem, output, surface_name = "Ground", decimate = 1L, boundary = NULL )export_landxml( dem, output, surface_name = "Ground", decimate = 1L, boundary = NULL )
dem |
A terra::SpatRaster representing the terrain surface. |
output |
Character. Output file path (should end in |
surface_name |
Character. Name for the TIN surface in the XML.
Default |
decimate |
Integer. Keep every Nth point to reduce file size.
Default |
boundary |
An optional sf::sf polygon to clip the DEM before export. |
The output file path (invisibly).
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) out <- file.path(tempdir(), "surface.xml") export_landxml(dem, out, surface_name = "Existing") file.exists(out)dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) out <- file.path(tempdir(), "surface.xml") export_landxml(dem, out, surface_name = "Existing") file.exists(out)
Converts a DEM raster to an STL (stereolithography) file suitable for 3D printing. The model includes a flat base and vertical exaggeration control.
export_stl(dem, output, exaggeration = 1, base_height = 1, decimate = 1L)export_stl(dem, output, exaggeration = 1, base_height = 1, decimate = 1L)
dem |
A terra::SpatRaster representing the terrain surface. |
output |
Character. Output file path (should end in |
exaggeration |
Numeric. Vertical exaggeration factor.
Default |
base_height |
Numeric. Thickness of the flat base below the
terrain minimum, in DEM units. Default |
decimate |
Integer. Keep every Nth cell. Default |
The output file path (invisibly).
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) out <- file.path(tempdir(), "terrain.stl") export_stl(dem, out, exaggeration = 2) file.exists(out)dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) out <- file.path(tempdir(), "terrain.stl") export_stl(dem, out, exaggeration = 2) file.exists(out)
Light-to-dark blue color ramp for flood depth visualisation.
flood_colors(n = 100L)flood_colors(n = 100L)
n |
Integer. Number of colors. Default |
Character vector of hex color values.
cols <- flood_colors(10) image(matrix(1:10), col = cols)cols <- flood_colors(10) image(matrix(1:10), col = cols)
Calculates the depth of water at each cell for a given water
surface elevation. Cells above the water level receive NA.
flood_depth(dem, water_level, boundary = NULL)flood_depth(dem, water_level, boundary = NULL)
dem |
A terra::SpatRaster representing the terrain surface. |
water_level |
Numeric. Water surface elevation. |
boundary |
An optional sf::sf polygon to restrict analysis to. |
A terra::SpatRaster of water depth values (positive,
in DEM elevation units). Cells above the water level are NA.
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) depth <- flood_depth(dem, water_level = 305) terra::plot(depth, main = "Flood Depth (m)")dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) depth <- flood_depth(dem, water_level = 305) terra::plot(depth, main = "Flood Depth (m)")
Creates a binary raster showing which cells would be inundated
at a given water surface elevation. Cells below the water level
are marked as flooded (1), cells above are NA.
flood_inundation(dem, water_level, boundary = NULL)flood_inundation(dem, water_level, boundary = NULL)
dem |
A terra::SpatRaster representing the terrain surface. |
water_level |
Numeric. Water surface elevation in the same units as the DEM. |
boundary |
An optional sf::sf polygon to restrict analysis to (e.g., a floodplain boundary). |
A terra::SpatRaster with values 1 (inundated) and
NA (dry).
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) inundated <- flood_inundation(dem, water_level = 305) terra::plot(inundated, main = "Flood Inundation at 305m")dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) inundated <- flood_inundation(dem, water_level = 305) terra::plot(inundated, main = "Flood Inundation at 305m")
Returns TRUE only when the ABOVER_KFA_TEST environment variable
is set, ensuring KyFromAbove examples never run on CRAN or in
environments without verified S3 access.
has_s3_access()has_s3_access()
Logical scalar indicating whether the KyFromAbove test environment variable is set.
has_s3_access()has_s3_access()
Calculates a relative elevation model representing each cell's height above the nearest drainage channel. This is a key input for flood susceptibility mapping. Channels are identified using Topographic Position Index (TPI).
height_above_drainage(dem, channel_threshold = 0.1, window = 15L)height_above_drainage(dem, channel_threshold = 0.1, window = 15L)
dem |
A terra::SpatRaster representing the terrain surface. |
channel_threshold |
Numeric. TPI percentile threshold for
channel identification. Lower values select deeper channels.
Default |
window |
Integer. Focal window size for TPI computation.
Must be odd. Default |
A terra::SpatRaster of height-above-nearest-drainage values. Lower values indicate greater flood susceptibility.
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) hand <- height_above_drainage(dem, window = 7) terra::plot(hand, main = "Height Above Nearest Drainage")dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) hand <- height_above_drainage(dem, window = 7) terra::plot(hand, main = "Height Above Nearest Drainage")
Computes an analytical hillshade (shaded relief) from a DEM using specified sun azimuth and altitude angles. Useful for terrain visualisation and map backgrounds.
hillshade(dem, azimuth = 315, altitude = 45)hillshade(dem, azimuth = 315, altitude = 45)
dem |
A terra::SpatRaster representing the terrain surface. |
azimuth |
Numeric. Sun azimuth in degrees (compass bearing).
Default |
altitude |
Numeric. Sun altitude angle in degrees above the
horizon. Default |
A terra::SpatRaster with hillshade values (0–255).
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) hs <- hillshade(dem) terra::plot(hs, col = grey.colors(256))dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) hs <- hillshade(dem) terra::plot(hs, col = grey.colors(256))
Calculates storage volume at a series of water surface elevations for a terrain depression (e.g., a pond, reservoir, or sediment basin). The result is an elevation-area-volume curve.
impoundment_curve(dem, boundary = NULL, elevations = NULL, n_steps = 20L)impoundment_curve(dem, boundary = NULL, elevations = NULL, n_steps = 20L)
dem |
A terra::SpatRaster representing the terrain surface. |
boundary |
An sf::sf polygon defining the impoundment boundary
(e.g., dam crest outline). If |
elevations |
Numeric vector of water surface elevations to
evaluate. If |
n_steps |
Integer. Number of elevation increments when
|
A data frame with columns:
elevation: water surface elevation
area_m2: inundated area at this elevation
volume_m3: cumulative storage volume below this elevation
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) boundary <- sf::st_read( system.file("extdata/boundary.gpkg", package = "aboveR"), quiet = TRUE ) curve <- impoundment_curve(dem, boundary, n_steps = 10) plot(curve$elevation, curve$volume_m3, type = "l", xlab = "Elevation", ylab = "Volume (m3)")dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) boundary <- sf::st_read( system.file("extdata/boundary.gpkg", package = "aboveR"), quiet = TRUE ) curve <- impoundment_curve(dem, boundary, n_steps = 10) plot(curve$elevation, curve$volume_m3, type = "l", xlab = "Elevation", ylab = "Volume (m3)")
Returns a bounding box (in EPSG:4326) for a Kentucky county by name. Useful for quickly querying KyFromAbove tiles without constructing an AOI manually.
kfa_county_bbox(county)kfa_county_bbox(county)
county |
Character. County name (case-insensitive). Partial matching is supported. |
A numeric vector c(xmin, ymin, xmax, ymax) in EPSG:4326.
# Get bounding box for Fayette County bbox <- kfa_county_bbox("Fayette") print(bbox) # Use directly with kfa_find_tiles tiles <- kfa_find_tiles(kfa_county_bbox("Fayette"), product = "dem")# Get bounding box for Fayette County bbox <- kfa_county_bbox("Fayette") print(bbox) # Use directly with kfa_find_tiles tiles <- kfa_find_tiles(kfa_county_bbox("Fayette"), product = "dem")
Queries KyFromAbove tile indexes to find elevation, point cloud, or imagery tiles that intersect a given area of interest. Tries the STAC catalog first (if available), then falls back to the tile index GeoPackage on S3.
kfa_find_tiles( aoi, product = "dem", phase = 2L, method = c("auto", "stac", "index") )kfa_find_tiles( aoi, product = "dem", phase = 2L, method = c("auto", "stac", "index") )
aoi |
An sf::sf object or numeric bbox ( |
product |
Character. One of |
phase |
Integer. KyFromAbove acquisition phase: 1, 2, or 3. Phase 1 DEMs are 5ft resolution; Phase 2/3 are 2ft. |
method |
Character. |
An sf::sf data frame with columns: tilename, s3_url, phase,
product, and geometry.
# Find Phase 2 DEM tiles for a bounding box in Fayette County tiles <- kfa_find_tiles( aoi = c(-84.55, 37.95, -84.45, 38.05), product = "dem", phase = 2 )# Find Phase 2 DEM tiles for a bounding box in Fayette County tiles <- kfa_find_tiles( aoi = c(-84.55, 37.95, -84.45, 38.05), product = "dem", phase = 2 )
Returns a character vector of all 120 Kentucky county names.
kfa_list_counties()kfa_list_counties()
Character vector of county names sorted alphabetically.
counties <- kfa_list_counties() head(counties)counties <- kfa_list_counties() head(counties)
Finds DEM tiles covering the AOI, reads them as Cloud-Optimized
GeoTIFFs via /vsicurl/, optionally merges into a single raster,
and crops to the AOI extent.
kfa_read_dem(aoi, phase = 2L, merge = TRUE, crop = TRUE, cache = FALSE)kfa_read_dem(aoi, phase = 2L, merge = TRUE, crop = TRUE, cache = FALSE)
aoi |
An sf::sf object or numeric bbox. |
phase |
Integer. 1 (5ft), 2 (2ft), or 3 (2ft). |
merge |
Logical. Mosaic multiple tiles into one SpatRaster? Default |
crop |
Logical. Crop result to AOI extent? Default |
cache |
Logical. Cache downloaded tiles locally? Default |
A terra::SpatRaster object.
dem <- kfa_read_dem( aoi = c(-84.55, 37.95, -84.45, 38.05), phase = 2 )dem <- kfa_read_dem( aoi = c(-84.55, 37.95, -84.45, 38.05), phase = 2 )
Finds ortho (nadir) or oblique imagery tiles covering the AOI and reads them as RGB SpatRaster.
kfa_read_ortho(aoi, phase = 3L, type = c("nadir", "oblique"))kfa_read_ortho(aoi, phase = 3L, type = c("nadir", "oblique"))
aoi |
An sf::sf object or numeric bbox. |
phase |
Integer. 1, 2, or 3. |
type |
Character. |
A terra::SpatRaster object with 3 bands (RGB).
ortho <- kfa_read_ortho( aoi = c(-84.55, 37.95, -84.54, 37.96), phase = 3 )ortho <- kfa_read_ortho( aoi = c(-84.55, 37.95, -84.54, 37.96), phase = 3 )
Finds point cloud tiles (LAZ for Phase 1, COPC for Phase 2/3)
covering the AOI and reads them via lidR::readLAS().
kfa_read_pointcloud(aoi, phase = 2L)kfa_read_pointcloud(aoi, phase = 2L)
aoi |
An sf::sf object or numeric bbox. |
phase |
Integer. 1, 2, or 3. |
A lidR::LAS object.
las <- kfa_read_pointcloud( aoi = c(-84.55, 37.95, -84.54, 37.96), phase = 2 )las <- kfa_read_pointcloud( aoi = c(-84.55, 37.95, -84.54, 37.96), phase = 2 )
Queries the KyFromAbove STAC catalog for items matching an area of interest and product type. Requires the rstac package.
kfa_stac_search(aoi, collection = NULL, datetime = NULL)kfa_stac_search(aoi, collection = NULL, datetime = NULL)
aoi |
An sf::sf object or numeric bbox. |
collection |
Character. STAC collection ID (e.g.,
|
datetime |
Character. ISO 8601 datetime or range
(e.g., |
The KyFromAbove STAC catalog is under active development by
Ian Horn (github.com/ianhorn/kyfromabove-stac). This function
will be activated when the catalog is live. In the meantime,
use kfa_find_tiles() with method = "index".
An sf::sf data frame of STAC items with asset URLs, or an error if the STAC catalog is not yet available.
# STAC catalog not yet available — use kfa_find_tiles() instead tiles <- kfa_find_tiles( aoi = c(-84.55, 37.95, -84.45, 38.05), product = "dem", phase = 2 )# STAC catalog not yet available — use kfa_find_tiles() instead tiles <- kfa_find_tiles( aoi = c(-84.55, 37.95, -84.45, 38.05), product = "dem", phase = 2 )
Downloads a tile index GeoPackage from the KyFromAbove S3 bucket
and caches it locally. Subsequent calls use the cached copy unless
it is older than max_age_days.
kfa_tile_index(product = "dem", phase = 2L, max_age_days = 30L)kfa_tile_index(product = "dem", phase = 2L, max_age_days = 30L)
product |
Character. One of |
phase |
Integer. KyFromAbove acquisition phase: 1, 2, or 3. |
max_age_days |
Integer. Re-download if cache is older than this. Default 30. |
An sf::sf data frame representing the tile index grid.
idx <- kfa_tile_index(product = "dem", phase = 2) head(idx)idx <- kfa_tile_index(product = "dem", phase = 2) head(idx)
Compares two DEMs of a pond or sediment basin to estimate the volume of accumulated sediment. The pond area is defined by a boundary polygon, and sedimentation is computed as fill (positive change) within that area.
pond_sedimentation(before, after, pond_boundary)pond_sedimentation(before, after, pond_boundary)
before |
A terra::SpatRaster of the pond before sedimentation. |
after |
A terra::SpatRaster of the pond after sedimentation. |
pond_boundary |
An sf::sf polygon defining the pond extent. |
A list with:
sediment_volume_m3: estimated sediment volume (positive fill)
mean_depth_change_m: mean elevation change within the pond
max_accumulation_m: maximum sedimentation depth
pond_area_m2: area of the pond boundary
change_raster: terra::SpatRaster of elevation change within the pond
before <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) after <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) boundary <- sf::st_read( system.file("extdata/boundary.gpkg", package = "aboveR"), quiet = TRUE ) sed <- pond_sedimentation(before, after, boundary) cat("Sediment volume:", sed$sediment_volume_m3, "m3\n")before <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) after <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) boundary <- sf::st_read( system.file("extdata/boundary.gpkg", package = "aboveR"), quiet = TRUE ) sed <- pond_sedimentation(before, after, boundary) cat("Sediment volume:", sed$sediment_volume_m3, "m3\n")
Compares a current DEM against a target (design grade) surface to quantify how much of a reclamation area has been restored to the desired elevation. Returns per-cell deviations and summary statistics.
reclamation_progress(current, target, boundary = NULL, tolerance = 0.3)reclamation_progress(current, target, boundary = NULL, tolerance = 0.3)
current |
A terra::SpatRaster of the current terrain surface. |
target |
A terra::SpatRaster of the target / design grade surface. |
boundary |
An optional sf::sf polygon to restrict analysis to. |
tolerance |
Numeric. Cells within +/- |
A list with:
deviation: terra::SpatRaster of current - target values
on_grade_pct: percentage of cells within tolerance
above_grade_pct: percentage of cells above tolerance
below_grade_pct: percentage of cells below tolerance
mean_deviation: mean signed deviation
rmse: root mean square error
current <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) target <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) prog <- reclamation_progress(current, target, tolerance = 1) cat("On grade:", prog$on_grade_pct, "%\n")current <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) target <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) prog <- reclamation_progress(current, target, tolerance = 1) cat("On grade:", prog$on_grade_pct, "%\n")
Calculates slope (in degrees or percent) and aspect (compass bearing) from a DEM in a single call.
slope_aspect(dem, units = c("degrees", "percent"))slope_aspect(dem, units = c("degrees", "percent"))
dem |
A terra::SpatRaster representing the terrain surface. |
units |
Character. Slope units: |
A two-layer terra::SpatRaster:
slope: terrain slope in the requested units
aspect: terrain aspect in degrees (0 = North, 90 = East, 180 = South,
270 = West; flat areas = -1)
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) sa <- slope_aspect(dem) terra::plot(sa)dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) sa <- slope_aspect(dem) terra::plot(sa)
Calculates surface roughness as the standard deviation of elevation within a moving window. Rougher surfaces indicate unreclaimed terrain, active construction, or natural heterogeneity.
surface_roughness(dem, window = 5L)surface_roughness(dem, window = 5L)
dem |
|
window |
Integer. Size of the moving window (number of cells).
Must be odd. Default |
A terra::SpatRaster of local roughness values (standard deviation of elevation).
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) rough <- surface_roughness(dem, window = 5) terra::plot(rough)dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) rough <- surface_roughness(dem, window = 5) terra::plot(rough)
Calculates the elevation difference between a before and after DEM, returning both the continuous change values and a classified layer (cut / stable / fill). Rasters are aligned automatically if their extents or resolutions differ (same CRS required).
terrain_change(before, after, tolerance = 0.1)terrain_change(before, after, tolerance = 0.1)
before |
A terra::SpatRaster representing the earlier DEM. |
after |
A terra::SpatRaster representing the later DEM. |
tolerance |
Numeric. Changes within +/- |
A two-layer terra::SpatRaster:
change: continuous elevation difference (after - before)
class: integer classification (1 = cut, 2 = stable, 3 = fill)
before <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) after <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) result <- terrain_change(before, after) terra::plot(result[["change"]])before <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) after <- terra::rast(system.file("extdata/dem_after.tif", package = "aboveR")) result <- terrain_change(before, after) terra::plot(result[["change"]])
Green-brown-white hypsometric color ramp for elevation maps.
terrain_colors(n = 100L)terrain_colors(n = 100L)
n |
Integer. Number of colors. Default |
Character vector of hex color values.
cols <- terrain_colors(10) image(matrix(1:10), col = cols)cols <- terrain_colors(10) image(matrix(1:10), col = cols)
Samples elevation values from a DEM at regular intervals along a transect line and returns a data frame of distance vs. elevation.
terrain_profile(dem, line, spacing = NULL)terrain_profile(dem, line, spacing = NULL)
dem |
A terra::SpatRaster representing the terrain surface. |
line |
An sf::sf object containing a single LINESTRING geometry defining the transect. Or a path to a GeoPackage/shapefile. |
spacing |
Numeric. Distance between sample points along the
line, in the CRS units of |
A data frame with columns:
distance: distance along the profile from the start point
elevation: sampled elevation value
x, y: coordinates of each sample point
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) line <- sf::st_read( system.file("extdata/profile_line.gpkg", package = "aboveR"), quiet = TRUE ) prof <- terrain_profile(dem, line) plot(prof$distance, prof$elevation, type = "l", xlab = "Distance", ylab = "Elevation")dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) line <- sf::st_read( system.file("extdata/profile_line.gpkg", package = "aboveR"), quiet = TRUE ) prof <- terrain_profile(dem, line) plot(prof$distance, prof$elevation, type = "l", xlab = "Distance", ylab = "Elevation")
Extracts summary statistics for each polygon zone, including min, max, mean, median, standard deviation, and cell count.
zonal_stats(r, zones, id_field)zonal_stats(r, zones, id_field)
r |
|
zones |
An sf::sf data frame of polygons defining analysis zones. |
id_field |
Character. Column name in |
An sf::sf data frame with columns:
zone identifier
min, max, mean, median, sd: summary statistics
cell_count: number of non-NA cells
area_m2: zone area
dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) zones <- sf::st_read( system.file("extdata/zones.gpkg", package = "aboveR"), quiet = TRUE ) zs <- zonal_stats(dem, zones, id_field = "zone_id") print(zs)dem <- terra::rast(system.file("extdata/dem_before.tif", package = "aboveR")) zones <- sf::st_read( system.file("extdata/zones.gpkg", package = "aboveR"), quiet = TRUE ) zs <- zonal_stats(dem, zones, id_field = "zone_id") print(zs)