#' PHI Calculation
#'
#' Calculate MBSS Physical Habitat Index (PHI), Paul et al. 2003.
#'
#' No special libraries are needed for this function.
#'
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Erik.Leppo@tetratech.com (EWL)
# 20170327
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#' @param fun.df Data frame of data to be processed. Must include the columns
#' specified in the parameter fun.names.
#' @param fun.names Column names; Optional. Defaults = c("SampID", "Strata",
#' "Area_acres", "Remote_020", "Shading_0100", "EpiSub_020", "BankStab_020",
#' "AvgRipWid_m50max", "InStrmHab_020", "InstrmWood_Num", "RiffQual_020",
#' "Embedd_0100", "Erosn_ExtR_075", "Erosn_ExtL_075", "Erosn_SevR_03",
#' "Erosn_SevL_03", "RoadDist_m")
# @param fun.SampID Unique identifier (Station/Site/Sample).
# @param fun.Strata Maryland Strata (Highland, EPiedmont, Coastal).
# @param fun.PHICalcYear 1994 or 2000 (2000 is the default but can ).
# @param fun.Area_acres Watershed area (acres).
# @param fun.Remote_020 Remoteness (0-20) [1994 only, Replaces with Road
# Distance].
# @param fun.Shading_0100 Percent shading (0-100).
# @param fun.EpiSub_020 Epifaunal Substrate (0-20).
# @param fun.BankStab_020 Bank Stability (0-20) [1994 only, replaced with
# Erosion].
# @param fun.AvgRipWid_m50max Average riparian zone width (meters
# , maximum of 50).
# @param fun.InStrmHab_020 Instream habitat (0-20).
# @param fun.InstrmWood_Num Tally of instream wood (number).
# @param fun.RiffQual_020 Riffle quality (0-20).
# @param fun.Embedd_0100 Percent embeddedness (0-100).
# @param fun.Erosn_ExtR_075 Erosion, extent, right bank (0-75) [2000 only,
# replaces bank stability].
# @param fun.Erosn_ExtL_075 Erosion, extent, left bank (0-75) [2000 only,
# replaces bank stability].
# @param fun.Erosn_SevR_03 Erosion, severity, right bank (0-3) [2000 only,
# replaces bank stability].
# @param fun.Erosn_SevL_03 Erosion, severity, left bank (0-3) [2000 only,
# replaces bank stability].
# @param fun.RoadDist_m Distance to nearest road (meters) [2000 only,
# replaces remoteness].
#' @return Returns a dataframe of the original data with the metric scoring and
#' PHI values.
#'
#' @examples
#' # calculate PHI
#' PHI <- PHIcalc(MBSS.PHI)
#' # examine data
#' #View(PHI)
#' # Example Save file
#' #datetime <- format(Sys.time(),"%Y%m%d_%H%M%S")
#' #write.table(PHI,paste0("PHI_",datetime,".tab"),row.names=FALSE,sep="\t")
#' @export
PHIcalc <- function(fun.df,fun.names=c("SampID"
, "Strata"
, "Area_acres"
, "Remote_020"
, "Shading_0100"
, "EpiSub_020"
, "BankStab_020"
, "AvgRipWid_m50max"
, "InStrmHab_020"
, "InstrmWood_Num"
, "RiffQual_020"
, "Embedd_0100"
, "Erosn_ExtR_075"
, "Erosn_ExtL_075"
, "Erosn_SevR_03"
, "Erosn_SevL_03"
, "RoadDist_m")){
# 0. QC ####
fun.df <- as.data.frame(fun.df)
# 0.1. QC column names
# Define dataframe names
names.default <- c("SampID"
, "Strata"
, "Area_acres"
, "Remote_020"
, "Shading_0100"
, "EpiSub_020"
, "BankStab_020"
, "AvgRipWid_m50max"
, "InStrmHab_020"
, "InstrmWood_Num"
, "RiffQual_020"
, "Embedd_0100"
, "Erosn_ExtR_075"
, "Erosn_ExtL_075"
, "Erosn_SevR_03"
, "Erosn_SevL_03"
, "RoadDist_m")
# decide if use input names or default.
# provided names in function call so if user does not provide the defaults are
# used.
msg <- paste0("\n Incorrect number of column names provided (n="
,length(fun.names),"). Default is 17.")
try(if(length(fun.names)!=length(names.default)) stop(msg))
# assign names
names.df <- fun.names
# create name parameters to be used in the code.
name.SampID <- names.df[1]
name.Strata <- names.df[2]
name.Area_acres <- names.df[3]
name.Remote_020 <- names.df[4]
name.Shading_0100 <- names.df[5]
name.EpiSub_020 <- names.df[6]
name.BankStab_020 <- names.df[7]
name.AvgRipWid_m50max <- names.df[8]
name.InStrmHab_020 <- names.df[9]
name.InstrmWood_Num <- names.df[10]
name.RiffQual_020 <- names.df[11]
name.Embedd_0100 <- names.df[12]
name.Erosn_ExtR_075 <- names.df[13]
name.Erosn_ExtL_075 <- names.df[14]
name.Erosn_SevR_03 <- names.df[15]
name.Erosn_SevL_03 <- names.df[16]
name.RoadDist_m <- names.df[17]
name.Erosn_SevR_02 <- "Erosn_SevR_02"
name.Erosn_SevL_02 <- "Erosn_SevL_02"
#
# 0.2. Error Checking field boundaries ####
# Not implemented. For now PHI will calculate with NA data not contributing.
# 20170428, ensure numeric values. if not get error
## will still get error if have a real "character" value (but that is ok)
fun.df[,names.df[3:17]] <- sapply(fun.df[,names.df[3:17]], as.numeric)
#
# 0.2.1. Strata (HIGHLAND, EPIEDMONT, COASTAL)
# Don't need to check. If not the proper names then no PHI will be calculated
# # 0.2.2. 0-3 parameters
# name.Erosn_SevR_03
# name.Erosn_SevL_03
# # 0.2.3. 0-20 parameters
# name.Remote_020
# name.EpiSub_020
# name.BankStab_020
# name.InStrmHab_020
# name.RiffQual_020
# # 0.2.4. 0-50 parameters
# name.AvgRipWid_m50max
# # 0.2.5. 0-75 parameters
# name.Erosn_ExtR_075
# name.Erosn_ExtL_075
# # 0.2.6. 0-100 parameters
# name.Shading_0100
# name.Embedd_0100
#
# 1.0. Data Manipulation.####
fun.df[, "Strata"] <- toupper(fun.df[, "Strata"])
# 1.1. Assign null (NA) values to extra columns
fun.df[,"Erosn_SevR_02"] <- NA
fun.df[,"Erosn_SevL_02"] <- NA
fun.df[,"BANKSTAB1994"] <- NA
fun.df[,"BANKSTAB2000"] <- NA
fun.df[,"REMOTE1994"] <- NA
fun.df[,"REMOTE2000"] <- NA
name.PHI <- c("REMOTE"
, "SHADING"
, "EPI"
,"BANKSTAB"
,"RIPWID"
,"INSTRHAB"
,"WOOD"
,"RIFFQUAL"
,"EMBEDD")
fun.df[,name.PHI] <- NA
fun.df[,"PHI.denom"] <- NA
# 1.2. Modify Values
# Replace Severity Values for use in calculations
fun.df[,name.Erosn_SevR_02] <- as.numeric(gsub(2,1.5
,fun.df[,name.Erosn_SevR_03]))
fun.df[,name.Erosn_SevR_02] <- as.numeric(gsub(3,2
,fun.df[,name.Erosn_SevR_02]))
fun.df[,name.Erosn_SevL_02] <- as.numeric(gsub(2,1.5
,fun.df[,name.Erosn_SevL_03]))
fun.df[,name.Erosn_SevL_02] <- as.numeric(gsub(3,2
,fun.df[,name.Erosn_SevL_02]))
#
# 2.0. Score Metrics (n=3, different for each Strata)####
# 2.1. Highlands
boo.Strata <- fun.df[,name.Strata]=="HIGHLAND"
fun.df[,"REMOTE1994"][boo.Strata==TRUE] <- fun.df[
,name.Remote_020][boo.Strata==TRUE]/20
fun.df[,"REMOTE2000"][boo.Strata==TRUE] <- (sqrt(fun.df[
,name.RoadDist_m][boo.Strata==TRUE])*0.733+0.615)/20
fun.df[,"SHADING"][boo.Strata==TRUE] <- (asin(sqrt(fun.df[
,name.Shading_0100][boo.Strata==TRUE]/100))-0.226)/1.171
fun.df[,"EPI"][boo.Strata==TRUE] <- (fun.df[
,name.EpiSub_020][boo.Strata==TRUE])/18
fun.df[,"BANKSTAB1994"][boo.Strata==TRUE] <- (sqrt(fun.df[
,name.BankStab_020][boo.Strata==TRUE]) - 1)/3.472
fun.df[,"BANKSTAB2000"][boo.Strata==TRUE] <- (sqrt(((fun.df[
,name.Erosn_ExtR_075][boo.Strata==TRUE]/(-15))*fun.df[
,name.Erosn_SevR_02][boo.Strata==TRUE])+
((fun.df[
,name.Erosn_ExtL_075][boo.Strata==TRUE]/(-15))*
fun.df[,name.Erosn_SevL_02][boo.Strata==TRUE]) + 20)-1)/3.472
fun.df[,"RIPWID"][boo.Strata==TRUE] <- (fun.df[
,name.AvgRipWid_m50max][boo.Strata==TRUE])/50
fun.df[,"INSTRHAB"][boo.Strata==TRUE] <- NA
fun.df[,"WOOD"][boo.Strata==TRUE] <- NA
fun.df[,"RIFFQUAL"][boo.Strata==TRUE] <- NA
fun.df[,"EMBEDD"][boo.Strata==TRUE] <- NA
fun.df[,"PHI.denom"][boo.Strata==TRUE] <- 5
#}
# 2.2. Piedmont
boo.Strata <- fun.df[,name.Strata]=="EPIEDMONT"
# # if(tolower(fun.df[,name.Strata])=="piedmont") {
fun.df[,"REMOTE1994"][boo.Strata==TRUE] <- fun.df[
,name.Remote_020][boo.Strata==TRUE]/16
fun.df[,"REMOTE2000"][boo.Strata==TRUE] <- (sqrt(fun.df[
,name.RoadDist_m][boo.Strata==TRUE])*0.733+0.615)/16
fun.df[,"SHADING"][boo.Strata==TRUE] <- ((asin(sqrt(fun.df[
,name.Shading_0100][boo.Strata==TRUE]/100))-
(1.7528-0.1990*(log10(fun.df[
,name.Area_acres][boo.Strata==TRUE]))))+1.142)/1.405
fun.df[,"EPI"][boo.Strata==TRUE] <- (fun.df[
,name.EpiSub_020][boo.Strata==TRUE] - 1)/17
fun.df[,"BANKSTAB1994"][boo.Strata==TRUE] <- (sqrt(fun.df[
,name.BankStab_020][boo.Strata==TRUE]) - 1)/3.243
fun.df[,"BANKSTAB2000"][boo.Strata==TRUE] <- (sqrt(((fun.df[
,name.Erosn_ExtR_075][boo.Strata==TRUE]/(-15))*fun.df[
,name.Erosn_SevR_02][boo.Strata==TRUE])+
((fun.df[,name.Erosn_ExtL_075][boo.Strata==TRUE]/(-15))*fun.df[
,name.Erosn_SevL_02][boo.Strata==TRUE]) + 20)-1)/3.243
fun.df[,"RIPWID"][boo.Strata==TRUE] <- NA
fun.df[,"INSTRHAB"][boo.Strata==TRUE] <- (fun.df[
,name.InStrmHab_020][boo.Strata==TRUE]-
(9.9876+1.5476*log10(fun.df[,name.Area_acres][boo.Strata==TRUE]))+
12.805)/15.745
fun.df[,"WOOD"][boo.Strata==TRUE] <- fun.df[
,name.InstrmWood_Num][boo.Strata==TRUE]/12
fun.df[,"RIFFQUAL"][boo.Strata==TRUE] <- (fun.df[
,name.RiffQual_020][boo.Strata==TRUE]-
(5.8467+2.4075*log10(fun.df[,name.Area_acres][boo.Strata==TRUE]))+
16.252)/19.637
fun.df[,"EMBEDD"][boo.Strata==TRUE] <- (100 - fun.df[
,name.Embedd_0100][boo.Strata==TRUE])/90
fun.df[,"PHI.denom"][boo.Strata==TRUE] <- 8
# # }
# 2.3. Coastal Plain
boo.Strata <- fun.df[,name.Strata]=="COASTAL"
# if(tolower(fun.df[,name.Strata])=="coastalplain") {
fun.df[,"REMOTE1994"][boo.Strata==TRUE] <- fun.df[
,name.Remote_020][boo.Strata==TRUE]/18.570
fun.df[,"REMOTE2000"][boo.Strata==TRUE] <- (sqrt(fun.df[
,name.RoadDist_m][boo.Strata==TRUE])*0.733+0.615)/18.570
fun.df[,"SHADING"][boo.Strata==TRUE] <- (asin(sqrt(fun.df[
,name.Shading_0100][boo.Strata==TRUE]/100))-0.226)/1.120
fun.df[,"EPI"][boo.Strata==TRUE] <- (fun.df[
,name.EpiSub_020][boo.Strata==TRUE]-
(3.5233+2.5821*log10(fun.df[,name.Area_acres][boo.Strata==TRUE]))+
13.199)/17.213
fun.df[,"BANKSTAB1994"][boo.Strata==TRUE] <- sqrt(fun.df[
,name.BankStab_020][boo.Strata==TRUE])/4.472
fun.df[,"BANKSTAB2000"][boo.Strata==TRUE] <- (sqrt(((fun.df[
,name.Erosn_ExtR_075][boo.Strata==TRUE]/(-15))*fun.df[
,name.Erosn_SevR_02][boo.Strata==TRUE])+((fun.df[
,name.Erosn_ExtL_075][boo.Strata==TRUE]/(-15))*fun.df[
,name.Erosn_SevL_02][boo.Strata==TRUE]) + 20))/4.472
fun.df[,"RIPWID"][boo.Strata==TRUE] <- NA
fun.df[,"INSTRHAB"][boo.Strata==TRUE] <- (fun.df[
,name.InStrmHab_020][boo.Strata==TRUE]-(0.5505+4.2475*log10(fun.df[
,name.Area_acres][boo.Strata==TRUE]))+15.094)/18.023
fun.df[,"WOOD"][boo.Strata==TRUE] <- (fun.df[
,name.InstrmWood_Num][boo.Strata==TRUE]-(-12.24+8.8120*log10(fun.df[
,name.Area_acres][boo.Strata==TRUE]))+28.903)/33.803
fun.df[,"RIFFQUAL"][boo.Strata==TRUE] <- NA
fun.df[,"EMBEDD"][boo.Strata==TRUE] <- NA
fun.df[,"PHI.denom"][boo.Strata==TRUE] <- 6
# }##IF.fun.Strata.END
#
# 3.0. More Data Manipulations (post scoring)####
# 3.1.0. 1994 vs. 2000 parameters. 2000 is primary and used if not NA.
# Else use 1994.
# 3.1.1. REMOTE
fun.df[,"REMOTE"] <- ifelse(is.na(fun.df[,"REMOTE2000"])
,fun.df[,"REMOTE1994"]
,fun.df[,"REMOTE2000"])
# 3.1.2. BANKSTAB
fun.df[,"BANKSTAB"] <- ifelse(is.na(fun.df[,"BANKSTAB2000"])
,fun.df[,"BANKSTAB1994"]
,fun.df[,"BANKSTAB2000"])
# 3.2. Rescale from 0:1 to 0:100 and set min and max of 0 and 100
fun.df[,name.PHI] <- apply(fun.df[,name.PHI],c(1,2), function(x) 100*max(0
,min(1,x)))
#
# 4.0. Final PHI Score. Sum Metrics and divide by the proper denominator####
fun.df[,"PHI.sum"] <- rowSums(fun.df[,name.PHI],na.rm=TRUE)
fun.df[,"PHI"] <- fun.df[,"PHI.sum"]/fun.df[,"PHI.denom"]
#
# 5.0. Return result dataframe####
return(fun.df)
#
}##fun.PHI.END
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.