In FloodHydrology/InundationHydRology: Inundation Hydrology

The goal of this analysis is to estimate contemporary storage capacity for the Palmer Lab study sites on the Eastern Shore. This is an intitial analysis, and it will likley evolve over time.

A few quick notes on data:

1) All data is housed in the Palmer Lab subdirectory. (Choptank/Nate/Storage_Capacity)

2) The 1m DEM was downloaded from imap.maryland.gov

3) Wetland polygons were digitized by hand and were only meant to be a rough estimate of wetland coverage. A more detailed analysis [that incoporates surveyed cross sections] will be conducted at a later date.

4) This script will be housed at floodhydrology.com

Step 1: Workspace Organization

As always, the first step is to define relevant directories, aquire required librairies, and download input data. Note, this anlaysis requires a DEM, wetland location, and wetland polygon, all of which can be found in the "Storage_Capacity" subdirectory on the Palmer Lab server. Further, the script also requires the user to specify the location of the "storage capacity" data directory, a scratch workspsace, and location of the Whitebox Tools executable. For SESYNC users, note that Whitebox cannot be initiated from the RStudio Server, and you will need to use the WynAnalytics virtual machine.

#Clear memory (KH says not to do this...maybe I'll convert eventually)
rm(list=ls(all=TRUE))

#Defin relevant working directories
data_dir<-"//storage.research.sesync.org/palmer-group-data/choptank/Nate/Storage_Capacity/"
scratch_dir<-"C:\\ScratchWorkspace\\"
wbt_dir<-"C:/WBT/whitebox_tools"

#Download packages 
library(raster)
library(sp)
library(rgdal)
library(rgeos)
library(dplyr)
library(parallel)
library(devtools)

#Download package from GIT
#install_github("FloodHydrology/InundationHydRology")
library(InundationHydRology)

#Download relevant data 
dem_jl<-raster(paste0(data_dir,"II_work/jl_dem"))
dem_jr<-raster(paste0(data_dir,"II_work/jr_dem"))
pnts<-readOGR(paste0(data_dir,"II_Work/."),"Wetland_Locations", verbose = F)
burn<-readOGR(paste0(data_dir,"II_Work/."),"wetland_burn")

#Add Unique ID to points
pnts$WetID<-seq(1,length(pnts))
burn$WetID<-sp::over(burn, pnts)$WetID

Step 2: Burn Wetlands into DEM

Burn wetlands into dem using the GIW_burn function. Note, here we are literally burning in a single elevation. [The wetland will look like a cylinder!!!] Eventually, we'll want to incorporate the surveyed cross sections to represent more complicated geometry.

#Process Jackson Ln property
mask_jl<-extent(dem_jl)
burn_jl<-crop(burn, mask_jl)
dem_burn_jl<-GIW_burn(dem_jl, burn_jl[1,])
for(i in 2:length(burn_jl)){
  dem_burn_jl<-GIW_burn(dem_burn_jl, burn_jl[i,])
}

#Process Jones Rd property
mask_jr<-extent(dem_jr)
burn_jr<-crop(burn, mask_jr)
dem_burn_jr<-GIW_burn(dem_jr, burn_jr[1,])
for(i in 2:length(burn_jr)){
  dem_burn_jr<-GIW_burn(dem_burn_jr, burn_jr[i,])
}

Here are a few examples of what the burned DEMs look like :)

par(mar=c(0,0,0,0))
plot(crop(dem_burn_jr,gBuffer(burn_jr[2,],width=100, byid = T)))
plot(burn_jr[2,], add=T)

Step 3: Delineate Wetlands

Here, we use WhiteBox Tools to identify depressions in the DEM, then we will delineate individual subsheds from those depressions.

3.1 Wetland Depression Identification

#Identify depressions in DEM
giws_jl<-GIW_identification(
            dem=dem_burn_jl,                    #DEM in question
            min_size=100,                       #Minimum depression size (in map units)
            workspace="C:\\ScratchWorkspace\\", #Scratch Workspace
            wbt_path="C:/WBT/whitebox_tools")   #WBG toolbox location

giws_jr<-GIW_identification(
            dem=dem_burn_jr,                    #DEM in question
            min_size=100,                       #Minimum depression size (in map units)
            workspace="C:\\ScratchWorkspace\\", #Scratch Workspace
            wbt_path="C:/WBT/whitebox_tools")   #WBG toolbox location

Below are the identified depressions for the Jackson Lane property.

par(mar=c(0,0,0,0))
plot(giws_jl)

3.2 Wetland Subshed Delineation

#Create rapper function for subshed delineation
fun<-function(WetID){
  #Define point of interest
  pnt<-pnts[pnts$WetID==WetID,]

  #Define dem and depressions of interest
  if(gIntersects(pnt, burn_jr)==T){
    dem=dem_burn_jr
    depressions=giws_jr
  }else{
    dem=dem_burn_jl
    depressions=giws_jl
  }

  #Execute subshed delineation function
  subshed<-GIW_subshed_delineation(
    dem = dem,
    depressions=depressions,
    wetland=pnt)

  #Write raster to folder
  writeRaster(subshed, 
              paste0(data_dir, "II_work/subshed_",pnt$WetID,".asc"),
              overwrite=T)
}

#Run function
lapply(seq(1,20), fun)

Step 4: Calculate Storage Capacity

Using method from Jones et al., 2018, we estimate contemporary storage capacity for each wetland.

4.1 Stage Storage Relationship

Create function to complete the storage capacity estimate below. Base it on the old dem, but include below surface inundation from PPR paper. then, print out each basin with inundation.

#Create rapper function for subshed delineation
fun<-function(WetID){
  #Define point of interest
  pnt<-pnts[pnts$WetID==WetID,]

  #Define dem and depressions of interest
  if(gIntersects(pnt, burn_jr)==T){
    dem=dem_burn_jr
  }else{
    dem=dem_burn_jl
  }

  #Define subshed 
  subshed<-raster(paste0(data_dir, "II_work/subshed_",WetID,".asc"))

  #Devleop stage storage relationship
  storage<-GIW_stage_storage(
          subshed = subshed, #Subshed delineated in previous step
          dem = dem,         #DEM
          z_max = 2,         #Maximum inundation depth
          dz = 1/12)         #Inundation Interval

  #Export
  storage$WetID<-WetID
  storage
}

#Run function
t0<-Sys.time()
cl <- makePSOCKcluster(detectCores()) #Create Clusters
clusterEvalQ(cl, library(raster))   
clusterEvalQ(cl, library(rgeos))
clusterExport(cl, c('pnts','burn_jr','dem_burn_jr',"dem_burn_jl", "data_dir",'GIW_stage_storage'), env=.GlobalEnv)  #Send Clusters function with the execute function
x<-parLapply(cl, seq(1,20), fun) #Run execute Function
stopCluster(cl)  #Turn clusters off
tf<-Sys.time()
tf-t0

#unlist output
storage<-do.call(rbind,x)

4.2 Spill point

Here, we estimate the spill point for each wetland.

#Create function to estimate the spill point and print plot 
fun<-function(WetID){
  #Define point of interest
  pnt<-pnts[pnts$WetID==WetID,]

  #Define dem and depressions of interest
  if(gIntersects(pnt, burn_jr)==T){
    dem=dem_burn_jr
  }else{
    dem=dem_burn_jl
  }

  #Define subshed 
  subshed<-raster(paste0(data_dir, "II_work/subshed_",WetID,".asc"))

  #Define Storage Curve
  storage<-storage[storage$WetID==WetID,]

  #Calculate max inundation
  inundation<-GIW_max_inundation(
    subshed,  #wateshed raster
    dem, #DEM for the analysis
    storage   #Storage Curve 
    )

  #Create layers for plotting
  subshed<-rasterToPolygons(subshed, dissolve=T)
  dem<-crop(dem, subshed)
  dem<-mask(dem, subshed)  

  #Estimate Below water surface (using W and Lane 2016)
  a_ha<-(storage$area[1]*(.3048^2))*0.0001
  dV_ham=0.25*(a_ha^1.4742)
  dV<-(dV_ham*10000)*(3.28084^3)
  storage$volume<-storage$volume+dV

  #Create Output Plots!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  png(paste0(scratch_dir,WetID,".png"), 
      width = 6, height=3, units="in", res=150)

  #Plot side by side
  par(mfrow=c(1,2))

  #Plot Inundation
  par(mar=c(1,1,1,1))
  plot(subshed, main=paste0("Wetland = ",pnt$Name))
  plot(dem, add=T, legend=F)
  plot(inundation, add=T, legend=F, col="blue")

  #Plot stage storage curve
  par(mar=c(4,4,1,1))
  par(mgp=c(2.1,1,0))
  plot(storage$z*12, storage$volume/gArea(subshed)*12, type="n",
       ps=12, cex.axis=10/12, cex.lab=14/12, 
       xlab="Relative Wetland Stage [in]", 
       ylab=expression("Specific Storage [in"^3*"/in"^2*"]"))
  abline(h=storage$volume[storage$outflow_length>0][1]/gArea(subshed)*12, 
         lty=2, lwd=2, col="grey30")
  points(storage$z*12, storage$volume/gArea(subshed)*12,
         pch=21, col="grey30", bg="grey70", cex=2)

  #Turn device off
  dev.off()

  #Create df of relevant storage capacity information 
  output<-data.frame(WetID=WetID,
                     area_subshed=gArea(subshed),
                     area_wetland=storage$area[storage$outflow_length>0][1], 
                     volume=storage$volume[storage$outflow_length>0][1])
  #Export!
  output
}

#Run function
output<-lapply(seq(1,20), fun)
output<-do.call(rbind,output)