Groundwater Contouring
Bill Rice
10/3/2019
R Markdown
The process of generating groundwater elevation contours assumes several things. The most important step in generating groundwater contours is developing a good conceptual understanding of the hydrogeology of a region under study. The key question is “Do all the well observations represent the same aquifer system or should they be considered seperate aquifer systems?” If you determine that the wells can be treated as a single aquifer system then a single groundwater contour map makes sense but if multiple aquifer systems are represented then it is most appropriate to generate a seperate groundwater contour map for each aquifer.
It is assumed that the user of this script has a strong conceptual understanding of hydrogeology.
This example assumes all wells are part of the same aquifer system.
Getting started with groundwater plots
The groundwater contouring script shown here includes many steps, functions, and several libraries.
First I will show the steps in creating groundwater contours using my data and at the end I will give example outputs from my data inputs.
Functions to make contour lines and spatial objects from well elevations
contour_function <- function(well_observations=wl_obs, year, season, crop_to_data = TRUE,datum="gs",map_output="return_raster"){
library("sp")
library("mapview")
library("raster")
library("fields")
library("maptools")
library("rgeos")
library("rgdal")
library("prettymapr")
library("raster")
library("rgeos")
library(lubridate)
# Subfunction 1: takes "season" input from outer function and converts to dates.
# Handles "fall", "spring", and specific date windows.
season_date_maker = function(year = 2015, season = "fall"){
if(season == "fall"){
fall_start = "09 01"; fall_end = "11 30" #format: "mm dd"
season_dates = c(ymd(paste(year, fall_start)), ymd(paste(year, fall_end)))
} else if (season == "spring"){
spring_start = "03 01"; spring_end = "05 31"
season_dates = c(ymd(paste(year, spring_start)), ymd(paste(year, spring_end)))
} else if(length(season) == 2 & is.Date(season)){season_dates = season
} else {print("Please specify fall, spring, or a 2-element array of date-type objects")}
return(season_dates)
}
gw_elevation_subsetting = function(year, season){
# dem_watershed <- raster(file.path(here::here(), "Reference_Spatial_Files_Temp" ,"merged_dem_siskiyou_10m_3310_brad.tif"))
# dem_watershed <- raster(file.path(here::here(), "Reference_Spatial_Files_Temp" ,"merged_dem_siskiyou_10m_3310_bill.tif"))
#Build season start and end dates from year/season or pull from specified dates
season_dates <- season_date_maker(year, season)
#### Subset water level data - first pass
gw_elevation.subset = subset(well_observations,
well_observations$date >= season_dates[1] &
well_observations$date <= season_dates[2] &
gs_to_ws_ft != "NA")
gw_elevation.subset <- merge(gw_elevation.subset, wells, by="well_code", no.dups=TRUE, all=FALSE)
# 2. Make point data spatial and reproject all spatial data ------------------
#### Make spatial point data
gw_elevation.spatial <- gw_elevation.subset
coordinates(gw_elevation.spatial) = ~longitude + latitude #Make the point-perforations table spatial
proj4string(gw_elevation.spatial) <- CRS("+init=epsg:4326") #assign WGS84 projection to coordinates
#Reproject wl points and basin boundary
gw_elevation.spatial <- spTransform(gw_elevation.spatial, crs(wsh))
# Check if the point is over the watershed
gw_elevation.spatial <- gw_elevation.spatial[wsh,]
gw_elevation.spatial$gp_elev_ft <- extract(dem_watershed, gw_elevation.spatial, method='bilinear')*3.28084
# Check how different the DEM is from the RP elevation
# gw_elevation.spatial$feetmsl_from_dem <- extract(dem_watershed, gw_elevation.spatial, method='bilinear')*3.28084
# gw_elevation.spatial$difference_vs_dem <- gw_elevation.spatial$feetmsl_from_dem - as.numeric(gw_elevation.spatial$gp_elev_ft)
# View(gw_elevation.spatial@data)
# hist(gw_elevation.spatial$difference_vs_dem, breaks=40)
# Remove elevations if RP and DEM are more than +/- 8 feet
# gw_elevation.spatial <- subset(gw_elevation.spatial, gw_elevation.spatial$difference_vs_dem <=8 &
# gw_elevation.spatial$difference_vs_dem >=-8)
#Subset data - set more specific date window
season_spread = 7 # 0 => 1-day season; 1 => 3-day season;
gw_elevation.spatial <- subset(gw_elevation.spatial,
gw_elevation.spatial$date >= (as.Date(median(gw_elevation.spatial$date))
- as.difftime(season_spread, unit="days")) &
gw_elevation.spatial$date <= (as.Date(median(gw_elevation.spatial$date))
+ as.difftime(season_spread, unit="days")) )
return(gw_elevation.spatial)
}
if("return_raster" %in% map_output){
gw_elevation.spatial <- gw_elevation_subsetting(year, season)
# Pre-crop the DEM
if(crop_to_data == TRUE){
gw_elevation.convex <- gConvexHull(gw_elevation.spatial)
gw_elevation.convex <- gUnaryUnion(gw_elevation.convex)
gw_elevation.convex <- gBuffer(gw_elevation.convex, width=500)
dem_contour_boundary <- crop(dem_watershed, gw_elevation.convex)
} else {dem_contour_boundary <- dem_watershed}
if(datum=="msl"){
#generate Thin Plate Spline model object
gw_elevation_subsetting(year, season)
gw_elevation.tps <- Tps(gw_elevation.spatial@coords, as.numeric(gw_elevation.spatial$wse_ft), verbose=FALSE)
}
if(datum=="gs"){
gw_elevation_subsetting(year, season)
gw_elevation.tps <- Tps(gw_elevation.spatial@coords, as.numeric(gw_elevation.spatial$gs_to_ws_ft), verbose=FALSE)
}
#Interpolate spline to all points in basin DEM raster
gw_elevation.interpolate <- interpolate(dem_contour_boundary, gw_elevation.tps)
# Mask interpolated raster to groundwater basin boundary
# gw_elevation.interpolate <- mask(gw_elevation.interpolate, boundary_basin)
# to do : maybe take out above statement, since we crop to basin below
gw_elevation.raster <- gw_elevation.interpolate
#Crop Raster to the convex hull (largest polygon covered, + a buffer) of our water measurement data
if(crop_to_data == TRUE){
gw_elevation.convex <- gConvexHull(gw_elevation.spatial)
buf1 <- gBuffer(gw_elevation.convex, width=500, byid=TRUE)
buf2 <- gUnaryUnion(buf1)
buf <- gBuffer(buf2, width=500)
gw_elevation.raster <- mask(gw_elevation.raster, buf)
}
#### crop map to basin
gw_elevation.raster <- mask(crop(gw_elevation.raster,basin), basin)
}
if(map_output=="return_spatial"){
gw_elevation.spatial <- gw_elevation_subsetting(year, season)
return(gw_elevation.spatial)
}
if(map_output=="return_raster"){
return(gw_elevation.raster)
}
}groundwater_contour_figure = function(gw_elevation_raster, gw_elevation_spatial, datum, basin_name){
con <- function(test, yes, no, filename="", ...) {
stopifnot(inherits(test, "Raster"))
if (!inherits(no, "Raster")) {
x <- reclassify(test, rbind(c(0,no), c(1,NA)))
} else {
x <- mask(no, test, maskvalue=TRUE)
}
if (!inherits(yes, "Raster")) {
yes <- reclassify(test, rbind(c(1,yes), c(0,NA)))
}
cover(x, yes, filename=filename, ...)
}
# tmap_options(max.raster = c(plot = 20000000, view = 20000000))
if(datum=="msl"){
#map colors to elevations (ftamsl) for each basin
if(basin_name == "Scott River Valley"){
cuts = c(2680, 2700, 2720, 2740, 2760, 2780, 2800, 2820, 2840, 2890, 2940)}
if(basin_name == "Shasta Valley"){cuts = seq(2450, 2850, 50)}
if(basin_name == "Butte Valley"){cuts = seq(4100, 4240, 20)}
}
if(datum=="gs"){
#map colors to elevations (ftamsl) for each basin
if(basin_name == "Scott River Valley"){cuts = seq(0, 140, 20)}
if(basin_name == "Shasta Valley"){cuts = seq(0, 140, 20)}
if(basin_name == "Butte Valley"){cuts = seq(0, 150, 20)}
}
if(datum=="msl"){contour_pal <- colorRampPalette(c("red","orange", "yellow","white","blue"))}
if(datum=="gs"){contour_pal <- colorRampPalette(c("blue","white", "yellow","orange","red"))}
gw_elevation_raster <- con(gw_elevation_raster > 0.0 , gw_elevation_raster , 0)
gw_elevation.contour <- rasterToContour(gw_elevation_raster, levels = cuts)
contour_pal(length(cuts)-1)
tm_shape(hill_basin) +
tm_raster(palette = dem_shade, n = 100, legend.show = FALSE) +
tm_shape(basin, name = "Groundwater Basin",is.master = T, unit=mapunits) + tm_borders(color_basin, lwd = 1, lty=2) +
tm_shape(road_main, name = road_main$FULLNAME, unit=mapunits) + tm_lines(color_main_road, lwd = 3) + tm_text("FULLNAME")+
tm_shape(cities) +
tm_polygons(color_cities, lwd=1, alpha = 0.2) +
tm_shape(gw_elevation_raster) + tm_raster(alpha = 0.5, palette = contour_pal(length(cuts)-1), legend.show = FALSE)+
tm_shape(gw_elevation.contour) +
tm_iso(remove.overlap = TRUE, along.lines = TRUE, overwrite.lines = TRUE) + tm_text("level", size = 1.5) +
tm_shape(gw_elevation_spatial) +
tm_dots(col="black", size=0.5, alpha=0.2 )+
tm_text("gs_to_ws_ft", size=0.8, alpha = 0.8, auto.placement = TRUE)+
# tm_text("well_code", size=0.8, auto.placement = TRUE) +
tm_shape(cities_centroid, name = "Towns", unit=mapunits) +
# tm_symbols(color_cities, border.lwd=NA, size = 0.5) +
tm_text("name", size = 0.8, xmod = cities_centroid$label_xmod, ymod = cities_centroid$label_ymod, col=color_cities_text) +
tm_scale_bar()+
tm_compass(position = c("left", "top"), type = "4star", size = 2)+
tm_add_legend(type="line", lwd = c(1,1), col = c(color_watershed, color_basin),
labels = c("Watershed", "Groundwater Basin"),
lty=c(1,2)) +
tm_add_legend(type="fill", alpha=0.2, col=c(color_cities), labels=c("Cities"))+
# tm_add_legend(type="fill", alpha=0.2, col=c(contour_pal(length(cuts))), labels=c(paste(cuts, "Feet Below Ground")), reverse=FALSE)+
tm_add_legend(type="symbol", col = color_cities, border.lwd = NA, size = 0.5, labels = "Towns")+
tm_add_legend(type = "line", lwd = c(3), col = c(color_main_road), labels = c("Highway 97")) +
tm_layout(legend.bg.color = "white", legend.frame = T)
}
gw_elevation.raster_spring_2018 <- contour_function(wl_obs,
year = 2018,
season = "spring",
map_output = c("return_raster"),
crop_to_data = TRUE,
datum="gs")
gw_elevation.spatial_spring_2018 <- contour_function(wl_obs,
year = 2018,
season = "spring",
map_output = c("return_spatial"),
crop_to_data = TRUE,
datum="gs")
gw_figure <- groundwater_contour_figure(gw_elevation_raster=gw_elevation.raster_spring_2018,
gw_elevation_spatial=gw_elevation.spatial_spring_2018,
datum="gs",
basin_name="Butte Valley")
gw_figure 
ttm()
gw_figureSource data
The data used in this plot comes from multiple sources but is available through the California Open Data Portal.
Example prints of data
head(wl_obs)## stn_id well_code wlm_id date rp_elev_ft gp_elev_ft
## 1 22363 414796N1223568W001 675900 2008-10-30 3344.6 3343.6
## 2 22363 414796N1223568W001 675899 2008-04-04 3344.6 3343.6
## 3 22363 414796N1223568W001 675897 2007-03-28 3344.6 3343.6
## 4 22363 414796N1223568W001 675896 2006-10-19 3344.6 3343.6
## 5 22363 414796N1223568W001 675895 2006-04-20 3344.6 3343.6
## 6 22363 414796N1223568W001 675894 2005-10-13 3344.6 3343.6
## rdng_ws reading_rp_ft wse_ft rpe_wse gs_to_ws_ft wlm_qa_desc wlm_desc
## 1 0 242.9 3101.7 242.9 241.9 <NA> Unknown
## 2 0 242.1 3102.5 242.1 241.1 <NA> Unknown
## 3 0 241.2 3103.4 241.2 240.2 <NA> Unknown
## 4 0 241.3 3103.3 241.3 240.3 <NA> Unknown
## 5 0 242.0 3102.6 242.0 241 <NA> Unknown
## 6 0 242.8 3101.8 242.8 241.8 <NA> Unknown
## wlm_acc_desc wlm_org_id wlm_org_name
## 1 Water level accuracy is unknown 1 Department of Water Resources
## 2 Water level accuracy is unknown 1 Department of Water Resources
## 3 Water level accuracy is unknown 1 Department of Water Resources
## 4 Water level accuracy is unknown 1 Department of Water Resources
## 5 Water level accuracy is unknown 1 Department of Water Resources
## 6 Water level accuracy is unknown 1 Department of Water Resources
## comments coop_agency_org_id coop_org_name
## 1 <NA> 1 Department of Water Resources
## 2 <NA> 1 Department of Water Resources
## 3 <NA> 1 Department of Water Resources
## 4 <NA> 1 Department of Water Resources
## 5 <NA> 1 Department of Water Resources
## 6 <NA> 1 Department of Water Resources
## data_table_source
## 1 DWR_Periodic_GWL
## 2 DWR_Periodic_GWL
## 3 DWR_Periodic_GWL
## 4 DWR_Periodic_GWL
## 5 DWR_Periodic_GWL
## 6 DWR_Periodic_GWLhead(wl_obs)## stn_id well_code wlm_id date rp_elev_ft gp_elev_ft
## 1 22363 414796N1223568W001 675900 2008-10-30 3344.6 3343.6
## 2 22363 414796N1223568W001 675899 2008-04-04 3344.6 3343.6
## 3 22363 414796N1223568W001 675897 2007-03-28 3344.6 3343.6
## 4 22363 414796N1223568W001 675896 2006-10-19 3344.6 3343.6
## 5 22363 414796N1223568W001 675895 2006-04-20 3344.6 3343.6
## 6 22363 414796N1223568W001 675894 2005-10-13 3344.6 3343.6
## rdng_ws reading_rp_ft wse_ft rpe_wse gs_to_ws_ft wlm_qa_desc wlm_desc
## 1 0 242.9 3101.7 242.9 241.9 <NA> Unknown
## 2 0 242.1 3102.5 242.1 241.1 <NA> Unknown
## 3 0 241.2 3103.4 241.2 240.2 <NA> Unknown
## 4 0 241.3 3103.3 241.3 240.3 <NA> Unknown
## 5 0 242.0 3102.6 242.0 241 <NA> Unknown
## 6 0 242.8 3101.8 242.8 241.8 <NA> Unknown
## wlm_acc_desc wlm_org_id wlm_org_name
## 1 Water level accuracy is unknown 1 Department of Water Resources
## 2 Water level accuracy is unknown 1 Department of Water Resources
## 3 Water level accuracy is unknown 1 Department of Water Resources
## 4 Water level accuracy is unknown 1 Department of Water Resources
## 5 Water level accuracy is unknown 1 Department of Water Resources
## 6 Water level accuracy is unknown 1 Department of Water Resources
## comments coop_agency_org_id coop_org_name
## 1 <NA> 1 Department of Water Resources
## 2 <NA> 1 Department of Water Resources
## 3 <NA> 1 Department of Water Resources
## 4 <NA> 1 Department of Water Resources
## 5 <NA> 1 Department of Water Resources
## 6 <NA> 1 Department of Water Resources
## data_table_source
## 1 DWR_Periodic_GWL
## 2 DWR_Periodic_GWL
## 3 DWR_Periodic_GWL
## 4 DWR_Periodic_GWL
## 5 DWR_Periodic_GWL
## 6 DWR_Periodic_GWLdem_watershed## class : RasterLayer
## dimensions : 2177, 1632, 3552864 (nrow, ncol, ncell)
## resolution : 23.2, 30.8 (x, y)
## extent : -182358, -144495.6, 382879.4, 449931 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : memory
## names : USGS_NED_one_arcsec_SiskiyouExtent3310
## values : 818.9161, 2829.046 (min, max)wsh## class : SpatialPolygonsDataFrame
## features : 1
## extent : -182349.4, -144492.3, 382882.8, 451771 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## variables : 15
## names : OBJECTID, TNMID, MetaSource, SourceData, SourceOrig, SourceFeat, LoadDate, GNIS_ID, AreaAcres, AreaSqKm, States, HUC8, Name, Shape_Leng, Shape_Area
## value : 1, {522489C1-3A24-4D0F-8F85-83A71EB7D499}, NA, NA, NA, NA, 2012/06/11, 0, 388548.36, 1572.4, CA,OR, 18010205, Butte, 2.33088722827, 0.170162999863basin## class : SpatialPolygonsDataFrame
## features : 1
## extent : -178047.4, -154063.6, 416200.7, 444544.4 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## variables : 7
## names : OBJECTID, Basin_Numb, Basin_Subb, Basin_Name, Basin_Su_1, Region_Off, GlobalID
## value : 399, 1-003, 1-003, BUTTE VALLEY, BUTTE VALLEY, NRO, {374607EF-33B2-4FA6-8118-C5C1C0B050D8}cities## class : SpatialPolygonsDataFrame
## features : 4
## extent : -167802, -158731.9, 396901.3, 441209 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## variables : 16
## names : STATEFP, PLACEFP, PLACENS, GEOID, NAME, NAMELSAD, LSAD, CLASSFP, PCICBSA, PCINECTA, MTFCC, FUNCSTAT, ALAND, AWATER, INTPTLAT, ...
## min values : 06, 19584, 02408156, 0619584, Dorris, Dorris city, 25, C1, N, N, G4110, A, 1819372, 0, +41.5789379, ...
## max values : 06, 78176, 02410347, 0678176, Tennant, Tennant CDP, 57, U1, N, N, G4210, S, 628135, 40653, +41.9649654, ...hill_basin## class : RasterLayer
## dimensions : 920, 1034, 951280 (nrow, ncol, ncell)
## resolution : 23.2, 30.8 (x, y)
## extent : -178042.8, -154054, 416205, 444541 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0.1856007, 0.9998591 (min, max)color_cities## [1] "#d7922d"cities_centroid## class : SpatialPointsDataFrame
## features : 4
## extent : -167197.3, -159499.4, 397686.6, 440097 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## variables : 3
## names : name, label_xmod, label_ymod
## min values : Dorris, -1.2, -0.5
## max values : Tennant, 2.2, -0.5