Commit 5d3e05c8 authored by Asquith, William H.'s avatar Asquith, William H.

initial commit

parent 3c62678d
# README (path ./visGWDBmrva/inst/miscR/README.md)
#### Author: William H. Asquith
#### Point of contact: William H. Asquith (wasquith@usgs.gov)
***
***
# DISCUSSION
This README file describes the contents of the above path. One or more **R** language scripts are provided in this directory. These scripts are not ever to be integral to the operation of the **visGWDB** software, but instead are present here as utilities or example parent scripts that call the **visGWDB** software. The author's general style is to continue to have the current working directory set to that holding `visGWDB.R` but with one of these scripts below open and able to be sourced from the **RStudio** console.
1. Script `mudge_in_streamgages_for_study_nearby_aquifer_levels.R` --- This script uses 16 streamgages through their identification numbers and treats in sequel each streamgage was if it were a well having all missing water-level measurements. But the feature here is that the neighborhood of groundwater levels can be accessed, hydrographs plotted, and the regional generalized additive model and support vector machines still used to produce monthly outputs. Thus, for an alluvial aquifer, one has the ability to isolate the groundwater database in space and time near a streamgage. A loop is involved and the `visGWDB.R` is called once for each streamgage. The streamgage is tacked onto the database, used as the `sites` variable for **visGWDB**, then removed from the database, and the next streamgage attended too.
library(dataRetrieval)
library(lubridate)
library(sp)
library(rgdal)
LATLONG <- CRS("+init=epsg:4269 +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")
# Albers Equal Area with settings
ALBEA <- CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
streamgages <- c("07287150", "07364133", "07364200", "07364150", "07369680",
"07370000", "07264000", "07288500", "07288280", "07367700",
"07368000", "07288650", "07077380", "07077555", "07047942",
"07369500")
streamgages <- sort(streamgages)
MUDGE <- dataRetrieval::readNWISsite(streamgages)
MUDGE[is.na(MUDGE$alt_va),]
MUDGE$alt_va[MUDGE$site_no == "07364133"] <- 150.54
#https://nationalmap.gov/epqs/pqs.php?x=-91.65611&y=33.86639&units=Feet&output=xml
#<USGS_Elevation_Point_Query_Service>
#<script/>
#<Elevation_Query x="-91.65611" y="33.86639">
#<Data_Source>3DEP 1/3 arc-second</Data_Source>
#<Elevation>150.54</Elevation>
#<Units>Feet</Units>
#</Elevation_Query>
#</USGS_Elevation_Point_Query_Service>
#
# But let us just use the elevation service for a land surface elevation
# throughout.
elevs <- c(
"07047942", "206.16",
"07077380", "250.39",
"07077555", "173.53",
"07264000", "222.71",
"07287150", "142.13",
"07288280", "140.03",
"07288500", "101.42",
"07288650", "92.44",
"07364133", "150.54",
"07364150", "128.58",
"07364200", "90.35",
"07367700", "77.15",
"07368000", "64.16",
"07369500", "78.41",
"07369680", "111.29",
"07370000", "56.99")
elevs <- as.numeric(elevs[seq(2,(length(elevs)),2)])
plot(MUDGE$alt_va, elevs)
MUDGE$ALT_VA <- elevs
MUDGE$SITE_BADGE <- MUDGE$site_no
MUDGE$DEC_LAT_VA <- MUDGE$dec_lat_va
MUDGE$DEC_LONG_VA <- MUDGE$dec_long_va
MUDGE$SPECIAL_REFALT_VA <- MUDGE$alt_va
xy <- sp::SpatialPoints(cbind(MUDGE$DEC_LONG_VA, MUDGE$DEC_LAT_VA))
slot(xy, "proj4string") <- LATLONG
xy_aea <- sp::spTransform(xy, ALBEA); en <- coordinates(xy_aea)
MUDGE$PROJ_EASTING <- en[,1]
MUDGE$PROJ_NORTHING <- en[,2]
if(! exists("GWmaster")) source("include/visGWDB_emplacer.R")
source("include/visGWDB_helpers.R") # to get the deparseDateTime
for(site in MUDGE$SITE_BADGE) {
message("ensuring ",site," is not in GWtrim")
GWtrim <- GWtrim[! GWtrim$SITE_BADGE == site, ]
}
myDT <- seq(as.POSIXct("1980-01-15", tz="UTC"),
as.POSIXct("2018-12-15", tz="UTC"), by="months")
n <- length(myDT)
myDATE <- deparseDateTime(myDT)
myYEAR <- deparseDateTimeYear(myDATE, want="year" )
myMONTH <- deparseDateTimeYear(myDATE, want="month" )
myDAY <- deparseDateTimeYear(myDATE, want="day" )
myDECADE <- deparseDateTimeYear(myDATE, want="decade")
myDAYS_IN_MONTH <- lubridate::days_in_month(myDT)
myFRAC_YEAR <- ((myMONTH-1) + myDAY/myDAYS_IN_MONTH)/12
# minus and plus this many feet
MUDGE$ylim.up <- MUDGE$ALT_VA + 50
MUDGE$ylim.dn <- MUDGE$ALT_VA - 150
USRCODE <- file("./user_code.R", "w")
cat("pozo.requested.prediction.dates <- seq(as.POSIXct('1980-01-15 12:00:00', tz='UTC'),\n", file=USRCODE)
cat(" as.POSIXct('2018-12-15 12:00:00', tz='UTC'),\n", file=USRCODE)
cat(" by='months')\n", file=USRCODE)
cat("tmp <- deparseDateTime(pozo.requested.prediction.dates)\n", file=USRCODE)
cat("pozo.year.limits <- range(deparseDateTimeYear(tmp, want='year'))\n", file=USRCODE)
cat("sites <- mudge.site\n", file=USRCODE)
cat("ylim.lock.to.this.interval <- mudge.ylim.lock.to.this.interval\n", file=USRCODE)
cat("postMloop <- FALSE\n", file=USRCODE)
close(USRCODE); rm(USRCODE)
for(site in MUDGE$SITE_BADGE) {
message(" Special site mudge: ", site)
GWtrim <- GWtrim[! GWtrim$SITE_BADGE == site, ] # insurance policy again of a leak
TMP_MUDGE <- MUDGE[MUDGE$SITE_BADGE == site,]
mudge.site <- site
mudge.ylim.lock.to.this.interval <- c(TMP_MUDGE$ylim.dn, TMP_MUDGE$ylim.up)
if(is.na(TMP_MUDGE$SPECIAL_REFALT_VA[1])) {
stop("critical failure, altitude (datum) of streamgage is NA")
}
Pozo <- GWtrim[1:n,]
Pozo[1:n,] <- rep(NA,length(Pozo[1,]))
Pozo$SITE_BADGE <- rep(site, n)
Pozo$DATE <- myDATE
Pozo$LEV_DT <- myDT
Pozo$YEAR <- myYEAR
Pozo$MONTH <- myMONTH
Pozo$DAY <- myDAY
Pozo$DECADE <- myDECADE
Pozo$DAYS_IN_MONTH <- myDAYS_IN_MONTH
Pozo$FRAC_YEAR <- myFRAC_YEAR
Pozo$DEC_LAT_VA <- rep(TMP_MUDGE$DEC_LAT_VA[1], n)
Pozo$DEC_LONG_VA <- rep(TMP_MUDGE$DEC_LONG_VA[1], n)
Pozo$ALT_VA <- rep(TMP_MUDGE$ALT_VA[1],n)
Pozo$LEV_VA <- rep( 0, n)
Pozo$LEV_AGE_CD <- rep("T", n)
Pozo$NAT_AQFR_CD <- rep("special_location", n)
Pozo$OTID_ID <- Pozo$NAT_AQFR_CD
Pozo$PROJ_EASTING <- rep(TMP_MUDGE$PROJ_EASTING[1], n)
Pozo$PROJ_NORTHING <- rep(TMP_MUDGE$PROJ_NORTHING[1], n)
GWtrim <- rbind(GWtrim, Pozo) # binding our "mudged" gage on to the bottom
# of the well database itself
source("visGWDB.R") # run the visGWDB software
GWtrim <- GWtrim[! GWtrim$SITE_BADGE == site, ] # yet another, insurance
# policy again of a leak, but sometimes useful to comment this line on in
# testing.
message("******************************************************************")
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment