| Title: | BC Aquifer Data Tools |
|---|---|
| Description: | Set of tools for processing BC Aquifer lithology and yield data. |
| Authors: | Steffi LaZerte [aut, cre] (ORCID: <https://orcid.org/0000-0002-7690-8360>), Christine Bieber [aut], David van Everdingen [aut], Andarge Baye [aut], Province of British Columbia [cph] |
| Maintainer: | Steffi LaZerte <[email protected]> |
| License: | Apache License (>= 2) |
| Version: | 1.0.0 |
| Built: | 2026-10-03 07:24:00 UTC |
| Source: | https://github.com/bcgov/bcaquiferdata |
This app allows you to load a shapefile and filter aquifer/well data according to region, explore data, and export cleaned files.
aq_app()aq_app()
aq_app()aq_app()
Removes data cache
cache_clean(bcmaps_cded = FALSE)cache_clean(bcmaps_cded = FALSE)
bcmaps_cded |
Logical. Whether or not to also remove CDED files cached with the bcmaps package. These are used by bcaquifertools for acquiring TRIM data, but may also be cached for use by other workflows. |
# cache_clean() # cache_clean(bcmaps_cded = TRUE)# cache_clean() # cache_clean(bcmaps_cded = TRUE)
This function downloads, updates or loads locally stored data. Currently this
function returns wells, wells_sf, or lithology data. Note that these
data are originally from GWELLS, but are cleaned and summarized for use in
the bcaquiferdata package. For example wells_sf is a spatial version of the
data, and lithology is a cleaned and standardized version of lithology.
wells also contains the new standardized lithology data, along with the
original lithology observations and intermediate classification steps to
simplify error tracing.
data_read(type, update = FALSE, permission = FALSE)data_read(type, update = FALSE, permission = FALSE)
type |
Character. Type of data to return, one of |
update |
Logical. Force update of the data? |
permission |
Logical. Permission to create the cache folder. If |
Under normal circumstances, users will not need to use this function as it is used internally by the main workflow functions. However, users may wish to overview entire datasets.
Bear in mind that the lithology cleaning and standardizing, while better than the original data, will almost certainly still have errors!
Data frame or spatial features object of the requested data.
wells <- data_read("wells")wells <- data_read("wells")
Update the GWELLs data stored locally.
data_update( type = c("wells", "lithology"), download = TRUE, permission = FALSE )data_update( type = c("wells", "lithology"), download = TRUE, permission = FALSE )
type |
Character. Type of data to update. One of "all", "wells", "lithology" |
download |
Logical. Whether to re-download and process the data
( |
permission |
Logical. Permission to create the cache folder. If |
data_update(type = "wells") data_update(type = "lithology")data_update(type = "wells") data_update(type = "lithology")
This function takes a shape file of a region and creates a DEM of the region.
Lidar data is stored locally as tiles. Tiles are only downloaded if they
don't already exist unless only_new = FALSE. TRIM data is obtained via the
bcmaps package and stored locally as tiles. Note: TRIM elevation is
coarser than Lidar. Use Lidar unless it is missing for your region of
interest.
dem_region( region, source = "lidar", out_file = NULL, overwrite = FALSE, buffer = 1, lidar_dir = NULL, only_new = TRUE, progress = httr::progress() )dem_region( region, source = "lidar", out_file = NULL, overwrite = FALSE, buffer = 1, lidar_dir = NULL, only_new = TRUE, progress = httr::progress() )
region |
sf simple features object. Shape file of the region of interest. |
source |
Character. Source of DEM, "lidar", "trim" or a file path (or vector of file paths) to a custom DEM file. See Details. |
out_file |
Character. File path of where to save tif of DEM clipped to
|
overwrite |
Logical. If |
buffer |
Numeric. Percent buffer to apply to the |
lidar_dir |
Character. File path of where Lidar tiles should be stored.
Defaults to the cache directory. Only applies when |
only_new |
Logical. Whether to download all Lidar tiles, or only new
tiles that don't exist locally. Defaults to TRUE. Only apples when |
progress |
Function. Progress bar to use. Generally leave as is. |
Lidar tiles are the newest tile available. If you have reason to need a
historical file, contact the bcaquiferdata team to discuss your use case.
stars spatiotemporal array object
Lidar data is obtained from the LidarBC portal. The tiles data frame
contains is an internally created data frame listing tiles and their
respective download locations. Tiles to download are selected based on
overlap between map tiles and the provided shapefile (region). These Lidar
tiles can be browsed and downloaded manually via the
LidarBC Open LiDAR Data Portal
The grid of map tiles is obtained from the BC Data Catalogue, BCGS 1:20,000 Grid
TRIM data is obtained via the bcmaps package from the BC government
Data Catalogue
based on overlap between map tiles and the provided shapefile (region).
If a file path or vector of file paths are provided, a local DEM is
loaded with stars and combined with stars::st_mosaic(). Note that it is
assumed the data is elevation in metres.
The fastest way to process these data is to use an stars proxy object which
keeps manipulations fast. However, to save a copy of the DEM, it is much
faster to use sf::gdal_utils() (rather than stars::write_stars() after).
Therefore, if out_file is NULL stars::st_mosaic() and sf::st_crop()
will be used to manipulate a stars proxy object. Alternatively, if out_file
is provided, sf::gdal_utils() will be used to combine, crop, and save a
copy of the Lidar or Trim DEM for the region provided. These means that
slightly different methods are used and it results in a very slightly
different cropped region. However, since the region file is buffered before
cropping in either method, it shouldn't affect well elevation calculations
downstream.
library(sf) library(stars) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Fetch Lidar DEM creek_lidar <- dem_region(creek_sf) plot(creek_lidar) # Fetch TRIM DEM creek_trim <- dem_region(creek_sf, source = "trim") plot(creek_trim) # Get dem and save to file creek_lidar2 <- dem_region(creek_sf, out_file = "clinton_lidar.tif") plot(creek_lidar) creek_trim <- dem_region(creek_sf, source = "trim", out_file = "clinton_trim.tif") plot(creek_trim) # Use local DEM koksilah_sf <- st_read("misc/data/Koksilah_watershed4/Koksilah_watershed4.shp") koksilah_dem <- dem_region( koksilah_sf, source = "misc/data/Koksilah_Watershed_DEM_2km_Buffer.tif")library(sf) library(stars) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Fetch Lidar DEM creek_lidar <- dem_region(creek_sf) plot(creek_lidar) # Fetch TRIM DEM creek_trim <- dem_region(creek_sf, source = "trim") plot(creek_trim) # Get dem and save to file creek_lidar2 <- dem_region(creek_sf, out_file = "clinton_lidar.tif") plot(creek_lidar) creek_trim <- dem_region(creek_sf, source = "trim", out_file = "clinton_trim.tif") plot(creek_trim) # Use local DEM koksilah_sf <- st_read("misc/data/Koksilah_watershed4/Koksilah_watershed4.shp") koksilah_dem <- dem_region( koksilah_sf, source = "misc/data/Koksilah_Watershed_DEM_2km_Buffer.tif")
This function is still in development. Use at your own risk!
Creates an Excel file with formulas for calculating drawdown. File is
pre-filled with wells within 100km of location (well or coordinates).
drawdown( location, rate = NA, duration = NA, overwrite = FALSE, file_name = NULL, update = FALSE )drawdown( location, rate = NA, duration = NA, overwrite = FALSE, file_name = NULL, update = FALSE )
location |
Numeric. Either a vector with longitude/latitude, or a well tag number. |
rate |
Numeric. The pumping rate in m3/day (defaults to |
duration |
Numeric. The duration of pumping in days (defaults to |
overwrite |
Logical. Overwrite existing file? |
file_name |
Character. Optional file name for created excel file. |
update |
Logical. Force update of the data? |
Creates excel file
drawdown(85199, rate = 343, duration = 180) drawdown(85199, rate = 343, duration = 180, overwrite = TRUE) drawdown(22966, rate = 343, duration = 180) drawdown(22966) drawdown(c(-123.5593, 48.647), rate = 3.97, duration = 180) drawdown(c(-123.5593, 48.647))drawdown(85199, rate = 343, duration = 180) drawdown(85199, rate = 343, duration = 180, overwrite = TRUE) drawdown(22966, rate = 343, duration = 180) drawdown(22966) drawdown(c(-123.5593, 48.647), rate = 3.97, duration = 180) drawdown(c(-123.5593, 48.647))
The flag_depth_mismatch flag identifies wells where the depth is not the
same as the depth of the final lithology layer. Fixing these well depths
means replacing the well depth with the depth of the final lithology layer.
fix_depth_mismatch(wells_sub)fix_depth_mismatch(wells_sub)
wells_sub |
sf spatial data frame. Subset of wells data output by
|
This is done automatically for leapfrog exports, but otherwise not.
Fixed wells_sub data frame
library(sf) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek_sf) # Fix well depths creek_wells_fixed <- fix_depth_mismatch(creek_wells) # Explore all fixes dplyr::select(creek_wells_fixed, dplyr::starts_with("fix"))library(sf) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek_sf) # Fix well depths creek_wells_fixed <- fix_depth_mismatch(creek_wells) # Explore all fixes dplyr::select(creek_wells_fixed, dplyr::starts_with("fix"))
The flag_yield_zero flag identifies wells where the yield was marked as 0
in GWELLS but probably should have been recorded as NA. Fixing these well
yields means replacing the yield of 0 with NA.
fix_yield_zero(wells_sub, fix = TRUE)fix_yield_zero(wells_sub, fix = TRUE)
wells_sub |
Data frame. The subsetted Wells data frame. |
fix |
Logical. Whether to apply the fix. |
A glossary of flag terms
flagsflags
flagsA data frame with 10 rows and 2 columns:
flag name
Flag description
Suggested solution to the problem
Clean and categorize lithology descriptions into primary, secondary, tertiary and final lithology categories. Generally this function is used internally when loading and cleaning GWELLS lithology.
lith_fix(desc = NULL)lith_fix(desc = NULL)
desc |
Character. Text string to convert/fix. |
However statements can be tested directly with this function to see how it works and for troubleshooting.
Data frame of lithology categorizations
lith_fix("sandy gravel") # basic spell checks lith_fix("saandy gravel")lith_fix("sandy gravel") # basic spell checks lith_fix("saandy gravel")
A spatial data frame of map tiles with corresponding links to Lidar tiles.
tilestiles
tilesA data frame with 7,129 rows and 5 columns:
Tile name
Spatial data
Projection
Lidar tile name
Link to Lidar tile
The spatial grid of map tiles is obtained from the BC Data Catalogue, BCGS 1:20,000 Grid
Links to Lidar tile urls are extracted from the list at the LidarBC Open LiDAR Data Portal
Example data of fixed wells subsetted and with elevation added
wells_eg_fixedwells_eg_fixed
wells_eg_fixedA data frame with 31 rows and 71 columns.
Example data of unfixed wells subsetted and with elevation added
wells_eg_unfixedwells_eg_unfixed
wells_eg_unfixedA data frame with 31 rows and 71 columns.
This function takes a region shape file and the DEM of a region (output of
dem_region()), subsets the wells data (from GWELLS) to this region and adds
the elevation data.
wells_elev(wells_sub, dem, dem_extra = NULL, update = FALSE)wells_elev(wells_sub, dem, dem_extra = NULL, update = FALSE)
wells_sub |
sf spatial data frame. Subset of wells data output by
|
dem |
stars simple features object. Output of |
dem_extra |
stars simple features object. Output of |
update |
Logical. Force update of the data? |
Because combining elevation data measured from different sources can
introduce artifacts, if two sources of elevation are provided (i.e. both
dem and dem_extra), the data will contain extra columns for assessment.
Specifically, there will be three elevation columns, rather than just one.
elev1 contains elevations from the primary source (dem), elev2 contains
elevations from the secondary source (dem_extra), elev contains the
combined elevation data, elev1 unless missing, then elev2.
sf spatial data frame
library(sf) library(ggplot2) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek_sf) # Fetch Lidar DEM creek_lidar <- dem_region(creek_sf) # Collect wells in this region with added elevation from Lidar creek_wells <- wells_elev(creek_wells, creek_lidar) ggplot() + geom_sf(data = creek_sf) + geom_sf(data = creek_wells, aes(colour = elev), size = 0.5, fill = "NA", show.legend = FALSE) + coord_sf(datum = st_crs(3005)) # BC Albers # OR Fetch TRIM DEM creek_trim <- dem_region(creek_sf, type = "trim") # Collect wells in this region with added elevation from Lidar creek_wells <- wells_elev(creek_wells, creek_trim) ggplot() + geom_sf(data = creek_sf) + geom_sf(data = creek_wells, aes(colour = elev), size = 0.5, fill = "NA", show.legend = FALSE) + coord_sf(datum = st_crs(3005)) # BC Albers # Dealing with missing data mill_sf <- st_read("misc/data/MillBayWatershed.shp") mill_wells <- wells_subset(mill_sf) mill_lidar <- dem_region(mill_sf) mill_trim <- dem_region(mill_sf, type = "trim") mill_wells <- wells_elev(mill_wells, dem = mill_lidar, dem_extra = mill_trim) # See how the elevation data is combined, `dem` (elev1) is the primary source. select(mill_wells, well_tag_number, elev1, elev2, elev) # Use local DEM koksilah_sf <- st_read("misc/data/Koksilah_watershed4/Koksilah_watershed4.shp") koksilah_wells <- wells_subset(koksilah_sf) koksilah_dem <- dem_region( koksilah_sf, source = "misc/data/Koksilah_Watershed_DEM_2km_Buffer.tif") koksilah_wells <- wells_elev(koksilah_wells, dem = koksilah_dem) # Plot p <- koksilah_wells |> st_transform(crs = st_crs(koksilah_dem)) plot(koksilah_dem, reset = FALSE, key.pos = NULL) plot(p["elev"], add = TRUE, pch = 20)library(sf) library(ggplot2) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek_sf) # Fetch Lidar DEM creek_lidar <- dem_region(creek_sf) # Collect wells in this region with added elevation from Lidar creek_wells <- wells_elev(creek_wells, creek_lidar) ggplot() + geom_sf(data = creek_sf) + geom_sf(data = creek_wells, aes(colour = elev), size = 0.5, fill = "NA", show.legend = FALSE) + coord_sf(datum = st_crs(3005)) # BC Albers # OR Fetch TRIM DEM creek_trim <- dem_region(creek_sf, type = "trim") # Collect wells in this region with added elevation from Lidar creek_wells <- wells_elev(creek_wells, creek_trim) ggplot() + geom_sf(data = creek_sf) + geom_sf(data = creek_wells, aes(colour = elev), size = 0.5, fill = "NA", show.legend = FALSE) + coord_sf(datum = st_crs(3005)) # BC Albers # Dealing with missing data mill_sf <- st_read("misc/data/MillBayWatershed.shp") mill_wells <- wells_subset(mill_sf) mill_lidar <- dem_region(mill_sf) mill_trim <- dem_region(mill_sf, type = "trim") mill_wells <- wells_elev(mill_wells, dem = mill_lidar, dem_extra = mill_trim) # See how the elevation data is combined, `dem` (elev1) is the primary source. select(mill_wells, well_tag_number, elev1, elev2, elev) # Use local DEM koksilah_sf <- st_read("misc/data/Koksilah_watershed4/Koksilah_watershed4.shp") koksilah_wells <- wells_subset(koksilah_sf) koksilah_dem <- dem_region( koksilah_sf, source = "misc/data/Koksilah_Watershed_DEM_2km_Buffer.tif") koksilah_wells <- wells_elev(koksilah_wells, dem = koksilah_dem) # Plot p <- koksilah_wells |> st_transform(crs = st_crs(koksilah_dem)) plot(koksilah_dem, reset = FALSE, key.pos = NULL) plot(p["elev"], add = TRUE, pch = 20)
Export wells data for use in Strater and Voxler
wells_export(wells_sub, id, type, dir = ".", zip = FALSE, preview = FALSE)wells_export(wells_sub, id, type, dir = ".", zip = FALSE, preview = FALSE)
wells_sub |
Data frame. Output of |
id |
Character. Id to prepend to all output files e.g., "id_lith.csv" |
type |
Character. Format in which to export. One of "strater", "voxler", "archydro", "leapfrog", or "surfer" (case-insensitive). |
dir |
Character. Directory where files should be exported to. Defaults to working directory. |
zip |
Logical. Whether to export a zip archive of the files to |
preview |
Logical. Whether to preview the exports ( |
If preview = FALSE, a vector of file names (or if zip = TRUE, a
single filename of the zipped archive); if preview = TRUE, a list of data
frames.
library(sf) # Load a shape file defining the region of interest creek <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek) # Fetch Lidar DEM creek_lidar <- dem_region(creek) # Collect wells in this region with added elevation from Lidar creek_wells <- wells_elev(creek_wells, creek_lidar) # Preview data for Strater p <- wells_export(creek_wells, id = "clinton", type = "strater", preview = TRUE) names(p) p[["strater_lith"]] p[["strater_collars"]] p[["strater_wells"]] # Export data for Strater wells_export(creek_wells, id = "clinton", type = "strater", zip = TRUE) # Export data for Voxler wells_export(creek_wells, id = "clinton", type = "voxler", zip = TRUE) # Export Arc Hydro wells_export(creek_wells, id = "clinton", type = "archydro", zip = TRUE) # Export Surver wells_export(creek_wells, id = "clinton", type = "surfer", zip = TRUE) wells_export(creek_wells, id = "clinton", type = "leapfrog", zip = TRUE)library(sf) # Load a shape file defining the region of interest creek <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek) # Fetch Lidar DEM creek_lidar <- dem_region(creek) # Collect wells in this region with added elevation from Lidar creek_wells <- wells_elev(creek_wells, creek_lidar) # Preview data for Strater p <- wells_export(creek_wells, id = "clinton", type = "strater", preview = TRUE) names(p) p[["strater_lith"]] p[["strater_collars"]] p[["strater_wells"]] # Export data for Strater wells_export(creek_wells, id = "clinton", type = "strater", zip = TRUE) # Export data for Voxler wells_export(creek_wells, id = "clinton", type = "voxler", zip = TRUE) # Export Arc Hydro wells_export(creek_wells, id = "clinton", type = "archydro", zip = TRUE) # Export Surver wells_export(creek_wells, id = "clinton", type = "surfer", zip = TRUE) wells_export(creek_wells, id = "clinton", type = "leapfrog", zip = TRUE)
Filter the GWELLS data returning only wells within the provided shapefile.
wells_subset( region, fix_bottom_intervals = TRUE, fix_depth_missing = TRUE, fix_yield_zero = TRUE, update = FALSE )wells_subset( region, fix_bottom_intervals = TRUE, fix_depth_missing = TRUE, fix_yield_zero = TRUE, update = FALSE )
region |
sf simple features object. Shape file of the region of interest. |
fix_bottom_intervals |
Logical. Whether to add 1m to bottom lithology
intervals that has no thickness (identified by |
fix_depth_missing |
Logical. Whether to fix missing well depths by
making them equal to the depth of the final lithology layer. Default
|
fix_yield_zero |
Logical. Whether to fix well yields of 0 by making them
|
update |
Logical. Force update of the data? |
library(sf) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek_sf)library(sf) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek_sf)
Yield records are extracted from lithology observations and added to the wells data.
wells_yield(wells_sub)wells_yield(wells_sub)
wells_sub |
sf spatial data frame. Subset of wells data output by
|
Data frame or sf spatial data frame with wells data and added yield from lithology.
library(sf) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek_sf) # Get yield data for these wells creek_yield <- wells_yield(creek_wells)library(sf) # Load a shape file defining the region of interest creek_sf <- st_read("misc/data/Clinton_Creek.shp") # Get wells within this region creek_wells <- wells_subset(creek_sf) # Get yield data for these wells creek_yield <- wells_yield(creek_wells)