Aquifer Thickness determined from ~300,000 well logs

The California DWR Online Well Completion Report Database contains ~1,000,000 digitized well logs. About 1/3 of those wells have information on the top and bottom of the perforated interval. Visualizing the mean and median depth of this parameter gives a sense of the aquifer in that region.

Rather than visualize 300,000 wells, observations were grouped into the hydrologic regions defined by DWR’s Bulletin 118.

library(sp)
library(XML)
library(rgdal)
library(rgeos)
library(raster)
library(scales)
library(tidyverse)
library(sf)
library(here)
library(maptools)
library(viridis)
library(leaflet)

s <- shapefile(here('shapefiles','I08_B118_CA_GroundwaterBasins.shp')) # Bulletin 118 GW subbasins
dat <- read_csv("OSWCR_201801.csv")
head(dat)

## # A tibble: 6 x 46
##   WCRNumber  LegacyLogNumber RegionOffice    CountyName LocalPermitAgency 
##   <chr>      <chr>           <chr>           <chr>      <chr>             
## 1 WCR2016-0… <NA>            DWR Northern R… Colusa     Colusa County Env…
## 2 WCR2016-0… e0306175        DWR Northern R… Colusa     Colusa County Env…
## 3 WCR2016-0… e0303164        DWR Northern R… Colusa     <NA>              
## 4 WCR2016-0… e0303934        DWR Northern R… Colusa     Colusa County Env…
## 5 WCR2016-0… e0306131        DWR Northern R… Colusa     Colusa County Env…
## 6 WCR2016-0… e0299482        DWR Northern R… Colusa     Colusa County Env…
## # ... with 41 more variables: PermitDate <chr>, PermitNumber <chr>,
## #   OwnerAssignedWellNumber <chr>, WellLocation <chr>, City <chr>,
## #   PlannedUseFormerUse <chr>, DrillerName <chr>,
## #   DrillerLicenseNumber <chr>, RecordType <chr>, DecimalLatitude <dbl>,
## #   DecimalLongitude <dbl>, MethodofDeterminationLL <chr>,
## #   LLAccuracy <chr>, HorizontalDatum <chr>, GroundSurfaceElevation <chr>,
## #   ElevationAccuracy <chr>, ElevationDeterminationMethod <chr>,
## #   VerticalDatum <chr>, Township <chr>, Range <chr>, Section <chr>,
## #   BaselineMeridian <chr>, APN <chr>, DateWorkEnded <chr>,
## #   WorkflowStatus <chr>, ReceivedDate <chr>, TotalDrillDepth <int>,
## #   TotalCompletedDepth <dbl>, TopOfPerforatedInterval <dbl>,
## #   BottomofPerforatedInterval <dbl>, CasingDiameter <dbl>,
## #   DrillingMethod <chr>, Fluid <chr>, StaticWaterLevel <dbl>,
## #   TotalDrawDown <int>, TestType <chr>, PumpTestLength <int>,
## #   WellYield <dbl>, WellYieldUnitofMeasure <chr>,
## #   OtherObservations <chr>, MTRS <chr>

Clean

• lat and long NAs
• lat/long that are waaaay out of line with reality

clean_dat <- dat %>%
  filter(!is.na(DecimalLatitude & !is.na(DecimalLongitude))) %>% # remove NAs
  filter(DecimalLatitude < 300 & DecimalLongitude > -1000) %>%
  mutate(long_fixed = ifelse(DecimalLongitude > 0, DecimalLongitude * -1, DecimalLongitude))

full <- clean_dat %>% filter(!is.na(TopOfPerforatedInterval) & !is.na(BottomofPerforatedInterval))

How many observations did we omit?

nrow(dat) - nrow(full)

## [1] 590231

Assign better names and calculate saturated thickness and extract only relevant columns for plotting.

full %>%
  mutate(b = BottomofPerforatedInterval - TopOfPerforatedInterval) %>%
  dplyr::select(long_fixed, DecimalLatitude, TopOfPerforatedInterval,
                BottomofPerforatedInterval, b) %>%
  rename(bot = BottomofPerforatedInterval, top = TopOfPerforatedInterval,
         lon = long_fixed, lat = DecimalLatitude) -> full_b

full_b %>% filter(bot < 2500 & top < 2500 & b > 0) -> full_b_reasonable

Convert spatial polygons and pts to sf objects.

gwb <- st_as_sf(s,
       coords = c("lon", "lat"), # for point data
       remove = F, # don't remove these lat/lon cols from df
       crs = 4326)

pts <- st_as_sf(full_b_reasonable,
                coords = c("lon","lat"),
                remove = F,
                crs = 4326)

In ArcMap, I:

1. imported DWR Bulletin 118 gw basins
2. imported pts as XY data in a GEOGRAPHIC reference system
3. projected XY data into WGS84 (same projection as groundwater basins)
4. Analysis Tools -> Intersection . Input pts and basins
5. Open attribute table, and export to .txt.
6. Bring into R

gwb_pts <- read_csv(here("shapefiles","gwb_pts_intersection.txt"))

Calculate

Group by basin and subbasin (Basin_Subb), calculate mean, median. Sd is useless because distributions aren’t normal.

gwb_pts %>% group_by(Basin_Subb) %>% summarise(median = median(b, na.rm = TRUE)) -> gwb_med
gwb_pts %>% group_by(Basin_Subb) %>% summarise(mean = mean(b, na.rm = TRUE)) -> gwb_mean

# join back into polygons df
gwb_full <- gwb %>%
  left_join(gwb_med, by="Basin_Subb") %>%
  left_join(gwb_mean, by = "Basin_Subb")

Map

bins <- c(0,50, 100,200,300,400,500,600)

pal <- colorBin("inferno", gwb_full$median, bins = bins) # viridis, magma, inferno, plasma

# put into right crs
gwb_full %>% st_transform(4326) -> gwb_full
gwb_full$Basin_Subb_2 <- gwb_full$Basin_Subb

# rename groups for legend
gwb_full %>% rename(median_perforated_interval_thickness = Basin_Subb,
                    mean_perforated_interval_thickness = Basin_Subb_2) -> temp

leaflet(width = "100%") %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = gwb_full, stroke = FALSE, smoothFactor = 0.2,
              color = ~pal(mean), label = ~paste("mean =", as.character(round(mean, 2))),
              fillOpacity = 0.8,
              group = "mean_perforated_interval_thickness") %>%
  addPolygons(data = gwb_full, stroke = FALSE, smoothFactor = 0.2,
              color = ~pal(median), label = ~paste("thickness = ", as.character(median),"ft."),
              fillOpacity = 0.8,
              group = "median_perforated_interval_thickness") %>%
  addLegend(pal = pal, values = gwb_full$median,
            title = ("thickness (ft.)")) %>%
  addLayersControl(overlayGroups = c("median_perforated_interval_thickness","mean_perforated_interval_thickness")) %>%
  hideGroup("median_perforated_interval_thickness")