################################################################################
# Script to import grain size distribution (GSD) date of bed material from NWIS
# and calculate median (d50) and 84th percentile (d84) grain sizes in [mm].
# 
# Need to create a "GSD" folder in the directory you are working within.
#
# One can check for the existence of this data at a list of sites here:
# http://nwis.waterdata.usgs.gov/nwis/qwdata?search_criteria=multiple_site_no&search_criteria=site_no_file_attachment&search_criteria=obs_count_nu&search_criteria=multiple_parameter_cds&submitted_form=introduction
#
# NOTE: code can be adjusted to download and process any NWIS WQ parameter.
#       Simply change the look up table to associate the parameter code with
#       the water quality constituent(s) you are interested in (if working with
#       multiple). This code can be simplified if you are working with only one
#       WQ parameter code (do not need a lookup table).
#
# Joel S Sholtes - PhD Candidate
# Colorado State University - Dept. Civil & Env. Engr
# May, 2014
# jsholtes@gmail.com
###############################################################################

## USGS Dry sieved bulk bed sediment sample parameter codes -> MOST COMMONLY AVAILABLE
## NOTE:  These GSD values will be finer than pebble counts b/c they are based on
##        Bulk sediment grabs.
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#  80164  - Bed sediment, sieve diameter, percent smaller than 0.0625 millimeters
#   to    - by = 2^x where x is an integer and = (-4:7) (phi scale)
#  80175  - Bed sediment, sieve diameter, percent smaller than 128 millimeters
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## USGS Direct measurement (pebble count?) bed sediment sample parameter codes
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#  69080  - Bed sediment, direct measurement, percent smaller than 2048 millimeters
#   to    - by = 2^(11:-4) where "by" doesn't seem to follow geom series
#  69140  - Bed sediment, direct measurement, percent smaller than 0.063 millimeters
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## USGS wet Seived bulk bed sediment sample
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 69600   - Sediment	Bed sediment, wet sieved, sieve diameter, percent smaller than 0.0625 millimeters
#  to     - by = 2^x where x is an integer and = (-4:7) (phi scale)
# 69612	  - Sediment	Bed sediment, wet sieved, sieve diameter, percent smaller than 128 millimeters
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## USGS Dry sieved bulk bed sediment sample (different size classes than series 1)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#  69032  - Bed sediment, direct measurement, percent smaller than 2048 millimeters
#   to    - by = 2^(11:-4) where "by" doesn't seem to follow geom series
#  69076  - Bed sediment, direct measurement, percent smaller than 0.063 millimeters
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

####################################################################################
## 1: import list of sites to download data for
####################################################################################

# Clear data
rm(list = ls())
# set working dir
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/SuspLoad")

# import table of USGS site numbers
site.table <- read.table("susp_sites_v10.csv", header = TRUE, sep = ",", 
                         stringsAsFactors = F)
sites <- site.table[c("site.no", "site.nm")]
# append GSD variables as columns to site table
sites[c("n", "d50.tot", "d50.sand","d50.min", "d50.max", "d50.med", "d50.avg", 
        "d85.min", "d85.max", "d85.med", "d85.avg")] <- NA

####################################################################################
## 2: Lookup table for grain size based on param code (choose one series)
####################################################################################

# ~~~~~~~~~~~~~~~~~~~~~~~(1) USGS Dry sieved bulk bed sediment~~~~~~~~~~~~~~~~~~~~~~~~~~
param.cd <- (80164:80175)
sed.size <- 2^(-4:7)
param.table <- data.frame(param = paste("p", param.cd, sep = ""),
                          diam = sed.size)
params <- "&parameter_cd=80164" # initialize the URL web query string

# ~~~~~~~~~~~~~~~~~~~~~~~(2) USGS Direct measurement (pebble count?)~~~~~~~~~~~~~~~~~~~~
param.cd <- (69080:69140)
sed.size <- c(2048,1536,1024,768,512,384,256,192,128,125,106,100,96,90,75,64,63,53,50,45,
           37.5,31.5,26.5,25,22.4,19,16,13.2,12.5,11.2,9.5,8,7.5,6.3,5.6,4.75,4,3.35,
           2.8,2.36,2,1.7,1.4,1.18,1,0.85,0.71,0.6,0.5,0.425,0.355,0.3,0.25,0.212,0.18,
           0.15,0.125,0.106,0.09,0.075,0.063)
param.table <- data.frame(param = paste("p", param.cd, sep = ""),
                          diam = sed.size)
params <- "&parameter_cd=69080" # initialize the URL web query string

# ~~~~~~~~~~~~~~~~~~~~~~~(3) USGS wet Seived bulk bed sediment sample~~~~~~~~~~~~~~~~~~~
param.cd <- (69600:69612)
sed.size <- c( 0.0625, 0.1, 2^(-3:7) )
param.table <- data.frame(param = paste("p", param.cd, sep = ""),
                          diam = sed.size)
params <- "&parameter_cd=69600" # initialize the URL web query string

# ~~~~~~~~~~~~~~~~~~~~~~~(4) USGS Dry sieved bulk bed sediment~~~~~~~~~~~~~~~~~~~~~~~~~~
param.cd <- (69032:69076)
sed.size <- c(2048,1536,1024,768,512,384,256,192,125,106,100,96,90,75,63,53,50,45,37.5,
              26.5,25,22.4,19,13.2,12.5,11.2,9.5,7.5,6.3,5.6,4.75,3.35,2.8,2.36,1.7,
              1.18,0.85,0.6,0.425,0.3,0.212,0.15,0.106,0.075,0.053)
param.table <- data.frame(param = paste("p", param.cd, sep = ""),
                          diam = sed.size)
params <- "&parameter_cd=69032" # initialize the URL web query string

####################################################################################
## 3: Once you have selected the appropriate GSD series, run script from here
####################################################################################

# create a concatenated string of web query param codes to paste in NWIS URL below
for (i in 2:length(param.cd)) {
  params0 <- paste("&parameter_cd=", param.cd[i], sep = "")
  params <- paste(params, params0, sep = "")
}

# This loop reads in a tab-delimited GSD file from NWIS webservice in [mm]
# then parses out the GSD data to approximate the d50 and d84 using
# linear interpolation between seive diameter classes
for (i in 1:nrow(sites)){
  site.no <- sites$site.no[i]
  # attach leading zero to site.no if necessary
  if(nchar(site.no) < 8){
    site.no <- paste("0", site.no, sep = "")
  }
  
  # some site IDs are made up (extended flow records) therefore no
  # no USGS peak data will be available
  if( !is.na( as.numeric(site.no) ) ) { 
  
  # creates web query URL string
  path <- paste("http://nwis.waterdata.usgs.gov/nwis/qwdata?qw_count_nu=1&data_type=qw&site_no=",
                site.no, "&begin_date=&end_date=&format=wide_rdb&agency_cd=USGS",
                params, sep = "")
  
  # Note that comment.char = "#" tells R to skip lines beginning with "#" as in
  # USGS RDB files.
  gsd <- read.table(path, comment.char="#", header = T, sep = "\t", 
                     stringsAsFactors = F, na.strings = "")
  gsd <- gsd[-1,] # removes first row, which does not contain data
  
  # Some sites may not have bed sediment samples; this skips them.
  if (nrow(gsd) > 0 & ncol(gsd) > 14) {
    
      # create numeric matrix of GSD data (from character values)
      gsd.mat <- data.matrix(gsd[,13:ncol(gsd)])
      
      # grain sizes start in 13th col in this case, may need to be changed for other
      # WQ parameters
      param.cds <- colnames(gsd.mat)
      
      # select grain sizes from 'param.table' where param codes exist in site data
      # this is like a 'lookup' table routine. Not every site has the entire range
      # of sieve diameters available
      diams <- param.table$diam[param.table$param %in% param.cds]
      
      # Create temporary table to hold d50 and d84 for each sample
      d.table <- data.frame(matrix(data = NA, nrow = nrow(gsd), ncol = 2),
                            row.names = NULL, stringsAsFactors = F)
      names(d.table) <-  c("d50", "d84")
  
      # this loop calcs d50 and d84 for each bed sediment sample using approx(),
      # which is linear interpolation. Returns 'NA' if not enough data available.
      # First, we check if any of the samples had a d50 > 0.0625mm (first col of gsd.mat)
      # if > 50% of sample is less than 0.0625 then we calc d50 based on sand material only.
      # in some cases this results in a D50 > sand size material, which is biased. 
      # Looking back at this, I don't agree with normalizing and just calcing d50 based on 
      # sand size and greater. I think the code should be changed to simply say d50 is 
      # less than the best available data or not calculate it at all.
      if (min(gsd.mat[,1], na.rm = T) > 50) {
        d50.sand <- T
        d50.tot <- "< 0.0625 mm"
        for (j in 1:nrow(gsd) ){
          gsd0 <- gsd.mat[j,which(!is.na(gsd.mat[j,]))]
          if(length(gsd0) - sum(is.na(gsd0)) > 2) {
            # ID value of percent < 0.0625 mm and normalize data to this to calc GSD of 
            # material >= 0.0625 mm.
            diams0 <- diams[which(!is.na(gsd.mat[j,]))] 
            p63 <- min(gsd0)
            gsd0 <- (gsd0 - p63) / max(gsd0 - p63) * 100
            
            d.table$d50[j] <- signif( approx(gsd0, diams0, xout = 50)$y, 3)
            d.table$d84[j] <- signif( approx(gsd0, diams0, xout = 84)$y, 3)
          }
        }
      } else {
        d50.sand <- F
        d50.tot <- "total"
          for (j in 1:nrow(gsd) ){
            gsd0 <- gsd.mat[j,which(!is.na(gsd.mat[j,]))]
            # are there more than 2 pieces of data for linear interpolation?
            if(length(gsd0) - sum(is.na(gsd0)) > 2) { 
              # Extract diam vals where numeric GSD vals exist
              diams0 <- diams[which(!is.na(gsd.mat[j,]))] 
              d.table$d50[j] <- signif( approx(gsd0, diams0, xout = 50)$y, 3)
              d.table$d84[j] <- signif( approx(gsd0, diams0, xout = 84)$y, 3)
            }
          }
        }
      
      # number of rows that contain data
      n <- sum( !is.na(d.table$d50) )
      
      # summary stats of d50 and d84
      d50.min <- min(d.table$d50, na.rm = T)
      d50.max <- max(d.table$d50, na.rm = T)
      d50.med <- median(d.table$d50, na.rm = T)
      d50.avg <- signif( mean(d.table$d50, na.rm = T), 3)
      
      d84.min <- min(d.table$d84, na.rm = T)
      d84.max <- max(d.table$d84, na.rm = T)
      d84.med <- median(d.table$d84, na.rm = T)
      d84.avg <- signif( mean(d.table$d84, na.rm = T), 3)
      
      # write summary stats to appended site table
      sites[i,3:13] <- c(n, d50.tot, d50.sand, d50.min, d50.max, d50.med, d50.avg, 
                         d84.min, d84.max, d84.med, d84.avg)
      
      # create data frame for GSD data and write table to "GSD" folder
      gsd.table <- data.frame(gsd.mat, row.names = NULL)
      names(gsd.table) <- diams
      write.csv(gsd.table, paste("GSD/" , site.no, ".csv", sep = ""), 
                row.names = F) # Need to make this folder first
    }
  }
}

# save summary stat results
write.csv(sites, "results_GSD_Susp.csv", row.names = F)
