Lake Volume and Area Calculation

Lake_VolArea.r

## Code was compiled by Paul Julian
## contact info: pauljulianphd@gmail.com

# GIS libraries 
library(sp)
library(rgdal)
library(rgeos)
library(raster)

# projection
utm17=CRS(SRS_string ="EPSG:26917")

# read geotiff into R
## R does not interface with geotiff geodatabases (ESRI proprietary extension)
## Must export geodatabase geotiff as standalone geotiff in ArcGIS

bath=raster::raster("...export_raster/Bathym_50ft.tif"))# input the path to your data
proj4string(bath)<-utm17 #set rasters projection
bath=bath-(-1.32) # converts from NAVD88 to NGVD29 based on SFWMD conversion layer
bath.m=bath*0.3048 # convert inches to meters

plot(bath.m) #explore

## Idea from 
## Jones CN, Evenson GR, McLaughlin DL, et al (2018) 
## Estimating restorable wetland water storage at landscape scales. 
## Hydrological Processes 32:305–313. 
## doi: https://doi.org/10.1002/hyp.11405

# Functions to calculate volume and area
Con<-function(condition,trueValue,falseValue){
  return(condition*trueValue+(!condition)*falseValue)
}

inundate<-function(z,temp.grd){
  area<-temp.grd<z # Con(temp.grd>z,0,1)
  volume<-((z-temp.grd)*area)*ft.to.m(raster::res(area)[1])*ft.to.m(raster::res(area)[2])
  rslt=data.frame(z=z,
                  area=cellStats(area,'sum')*ft.to.m(raster::res(area)[1])*ft.to.m(raster::res(area)[2]),
                  volume=cellStats(volume,'sum'))
  return(rslt)
}

ft.to.m<-function(x){x * 0.3048}

# Individual value
inundate(2.6,bath.m)

z.val=seq(8.8,18.6,0.2)
rslt=data.frame()
for(i in 1:length(z.val)){
tmp=inundate(ft.to.m(z.val[i]),bath.m)
rslt=rbind(tmp,rslt)
print(i)
}
rslt$area.km2=rslt$area*1e-6
rslt$volume.km3=rslt$volume*1e-9
rslt

I found an error in the area calculation. The gist is updated.