### Uses iNaturalist species data with <= 30 m positional accuracy for land cover
### utilization/availability analysis
### In order to reduce bias introduced by human population density and proximity to roadways,
### individual plant records were first weighted by a human population density index
### (within 20 km radius, 1 km resolution) and by the lowest value of the percent of random points
### near roads (< 20 m) among land cover classes. The National Land Cover Dataset (NLCD 2016)
### is sampled to obtain population adjusted frequency of land cover types for species
##
# Specify  directory for saved functions
#FunctDirect <- "D:/My Drive/win-library/"
FunctDirect <- "C:/Users/jamesltracy61/Documents/R/win-library/"
setwd(FunctDirect)
# Read in User Defined Functions
source("SpatialFilter_Function.R")
source("EvalStatVars_Summary_Function.R")
source("Round2_Function.R")
source("ClassContingencyTable2_Function.R")
source("Kappa.Calc_Function.R")  # used in ClassContingencyTable_Function.R
source("Round2_Function.R")
##
# Load needed packages
#install.packages(c("maptools", "rgdal", "sp", "raster", "plyr","foreign","dplyr","ggpubr","data.table","electivity","reshape","ggcorrplot","psych")) 
library(maptools)
library(rgdal)
library(sp)
library(raster)
library(plyr)
library(foreign)
library(dplyr)
library(ggpubr)
library(data.table)
library(electivity)
library(reshape2)
library(ggcorrplot)
library(psych)
###
InDirect <- "K:/Monarch/"
## Specify focal mean radius in km for population density raster in geographic wgs84 projection
#FocalDist <- 20
StudyAreaName <- paste0("SCStatesEcoregionDissGeo")
PopGridName <- paste0("sc2010pop1k")
setwd(InDirect)
Species <- paste0("A_cap")
SpecShort <- paste0("Acap")
LongSpecies <- "Asclepias_a_capricornu"
LongSpecies2 <- gsub("_", " ", LongSpecies)
iNatDate <- "1Jan2012_31Dec2020" # date of downloading iNaturalist data
TargetGroupList <- c("Asteraceae", "Apocynaceae", "Lamiaceae", "Verbenaceae")
# ("Solanaceae", "Boraginaceae", "Brassicaceae", "Caryophyllaceae", "Convolvulaceae", "Cucurbitaceae", "Fabaceae", "Geraniaceae", "Plantaginaceae", "Scrophulariaceae")
TargetGroupDates <- "1Jan2014_31Dec2018"
#
# Specify minimum percent land cover (without open water) allowed for land cover to be included in calculations
MinPerc <- 0.05   # Usually 0.01 (= 1%)
Year <- 2016  # year of NLCD data
SpatFiltBuff <- 1 # units in km for spatial filter
# Specify name of NLCD raster  in geographic WGS84 projection
NLCDRasterName <- "nlcd2016scstg"
NLCDMxtCoreTableName <- paste0(SpecShort, "100con1k_", Year, "nlcdstats.dbf")
# Specify study area shapefile for clipping data points
StudyAreaName <- "AcapMxtCoreSCPolyWGS84"
# Specify units of buffers
Units <- "km"
# Set projection
CRS.WGS84 <- CRS("+init=epsg:4326")
#
#################################################################################
## Read in raw iNaturalist data and reformat
setwd(InDirect)
SpeciesDataRaw.df <- data.frame(read.csv(paste0(LongSpecies, "_iNaturalist_", iNatDate, "Raw.csv"), stringsAsFactors=FALSE))
head(SpeciesDataRaw.df)
SpeciesDataRaw.df$Source <- "iNat"
#str(SpeciesDataRaw.df)
SpeciesDataRaw.df$Year <-  substr(SpeciesDataRaw.df$observed_on, nchar(SpeciesDataRaw.df$observed_on)-4+1, nchar(SpeciesDataRaw.df$observed_on))
SpeciesDataRaw.df$Year <- as.numeric(SpeciesDataRaw.df$Year)
#
SpeciesData.df1 <- SpeciesDataRaw.df[,c("Source", "id", "latitude", "longitude", "positional_accuracy", "user_id", "observed_on", "Year", "scientific_name", "place_guess", "url", "geoprivacy")]
head(SpeciesData.df1)
nrow(SpeciesData.df1)
# Keep only records with no geoprivacy
SpeciesData.df2 <- SpeciesData.df1[ which(SpeciesData.df1$geoprivacy==""),]
head(SpeciesData.df2)                                                                                                         
nrow(SpeciesData.df2)
# Delete records with any missing geocoordinates
SpeciesData.df2 <- SpeciesData.df2[complete.cases(SpeciesData.df2[,c("latitude","longitude")]),]
#str(SpeciesData.df2)
# Keep only records with positional accuracy of <= 30 m
SpeciesData.df3 <- SpeciesData.df2[ which(SpeciesData.df2$positional_accuracy<=30),]
head(SpeciesData.df3)
nrow(SpeciesData.df3)
# Keep only records within two years of NLCD Year
SpeciesData.df4 <- SpeciesData.df3[ which(SpeciesData.df3$Year>=(Year-2)&SpeciesData.df3$Year<=(Year+2)),]
head(SpeciesData.df4)
nrow(SpeciesData.df4)
#####################
# First, convert to spatial points data frame
xy <- SpeciesData.df4[,c("longitude","latitude")]
SpeciesData.spdf <- SpatialPointsDataFrame(coords = xy, data = SpeciesData.df4, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
head(SpeciesData.spdf)
nrow(SpeciesData.spdf)
##
#######################################
## Read in and process Target Group data for sample bias correction
##
GroupDataList <- list()
Count <- 0
for(Group in TargetGroupList) {
  #Group <- "Asteraceae"
  GroupDataRaw.df <- data.frame(read.csv(paste0(Group, "_iNaturalist_", TargetGroupDates, "Raw.csv"), stringsAsFactors=FALSE))
  head(GroupDataRaw.df)
  GroupDataRaw.df$Source <- "iNat"
  #str(GroupDataRaw.df)
  GroupDataRaw.df$Year <-  substr(GroupDataRaw.df$observed_on, nchar(GroupDataRaw.df$observed_on)-4+1, nchar(GroupDataRaw.df$observed_on))
  GroupDataRaw.df$Year <- as.numeric(GroupDataRaw.df$Year)
  #
  GroupDataRaw.df$family <- Group
  GroupData.df1 <- GroupDataRaw.df[,c("Source", "id", "latitude", "longitude", "positional_accuracy", "user_id", "observed_on", "Year", "scientific_name", "family", "place_guess", "url", "geoprivacy")]
  head(GroupData.df1)
  nrow(GroupData.df1)
  # Keep only records with no geoprivacy
  GroupData.df2 <- GroupData.df1[ which(GroupData.df1$geoprivacy==""),]
  head(GroupData.df2)
  nrow(GroupData.df2)
  # Delete records with any missing geocoordinates
  GroupData.df2 <- GroupData.df2[complete.cases(GroupData.df2[,c("latitude","longitude")]),]
  #str(GroupData.df2)
  # Keep only records with positional accuracy of <= 30 m
  GroupData.df3 <- GroupData.df2[ which(GroupData.df2$positional_accuracy<=30),]
  head(GroupData.df3)
  nrow(GroupData.df3)
  # Keep only records within two years of NLCD Year
  GroupData.df4 <- GroupData.df3[ which(GroupData.df3$Year>=(Year-2)&GroupData.df3$Year<=(Year+2)),]
  head(GroupData.df4)
  nrow(GroupData.df4)
  Count <- Count + 1
  GroupDataList[[Count]] <- GroupData.df4
}
GroupDataAll.df <- do.call(rbind, GroupDataList)
head(GroupDataAll.df)
tail(GroupDataAll.df)
nrow(GroupDataAll.df)
#
#####################
# First, convert to spatial points data frame
xy <- GroupDataAll.df[,c("longitude","latitude")]
GroupData.spdf <- SpatialPointsDataFrame(coords = xy, data = GroupDataAll.df, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
head(GroupData.spdf)
nrow(GroupData.spdf)
##
# Save unthinned presence points that are not yet clipped to study area
writeOGR(SpeciesData.spdf, ".", paste0(Species, "_HiRes_UnClipped_UnThinned"), driver="ESRI Shapefile", overwrite=TRUE)
writeOGR(GroupData.spdf, ".", paste0(Species, "_TargetGroupHiRes_UnClipped_UnThinned"), driver="ESRI Shapefile", overwrite=TRUE)

##########################################################################################################
## NOTE: Use Arc Python script "LandCoverSelectivityAssessment_PointandMaxEntModelNLCDProcessing_Acap.py"
## to create StudyAreaName shapefile, clip the species and target group shapefiles, spatially thin these
## files to 1 km, and add distance to roads. Also create 10,000 random
### points shapefile for 20km buffered minimum convex polygon of thinned presence points for study area and
### add distsance to roads. Finally, obtain percent area of each land cover type for
### each MaxEnt model clipped to the study area.
##
##########################################################################################################
##


##############################################################################
### In order to reduce sample bias introduced by human population density and proximity
### to roadways, individual plant records are first weighted by a human population
### density and then by distance to roads.
#############################################################################
##
setwd(InDirect)
## Read in population density grid
PopDenRaster <- raster(PopGridName)
######
## Read back in clipped and spatially thinned presence points
PresThinShp <- readOGR(".", paste0(Species, "_HiRes_Thinned", SpatFiltBuff,"km_rd"))  # NOTE: readShapePoly was giving error, so used readOGR instead
head(PresThinShp)
nrow(PresThinShp)
#
PresenceDat.spdf1 <- raster::extract(PopDenRaster, PresThinShp, sp=TRUE)
head(PresenceDat.spdf1)
# Convert back to data frame
PresenceDat.df <- as.data.frame(PresenceDat.spdf1, stringsAsFactors=FALSE)
head(PresenceDat.df)
# Delete last two columns
PresenceDat.df <- PresenceDat.df[,-c((ncol(PresenceDat.df)-1):(ncol(PresenceDat.df)))]
nrow(PresenceDat.df)
# Rename last column to "PopDen"
colnames(PresenceDat.df)[ncol(PresenceDat.df)] <- "PopDen"
names(PresenceDat.df)[names(PresenceDat.df) == "NEAR_DIST"] <- "RoadDist"
nrow(PresenceDat.df)
# Omit any rows with NA values for PopDen
# Make function to delete rows with NAs in a specified column
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}
PresenceDat.df <- completeFun(PresenceDat.df, "PopDen")
nrow(PresenceDat.df)
head(PresenceDat.df)
#
PresenceDat.df$RawCount <- 1
head(PresenceDat.df)
nrow(PresenceDat.df)
###################################################################################
### Add NLCD class names to data
#################################################################################
##
## Read in NLCD raster for desired region
NLCDRaster <- raster(NLCDRasterName)
#
## Read in National Land Cover Dataset Class Names
NLCDClassNames.df <- data.frame(read.csv("NLCDClassNames.csv"), stringsAsFactors=FALSE)
head(NLCDClassNames.df)
nrow(NLCDClassNames.df)
NLCDClassNames.df$ClassName <- as.character(NLCDClassNames.df$ClassName)
#
########################################################
## Read in 10,000 random points in  minimum convex polygon around presence points
RandomShp <- readOGR(".", paste0(Species, "_RandomPoints_Thinned1km_10000rd"))
head(RandomShp)
# Extract PopDen data
RandomDat.spdf1 <- raster::extract(PopDenRaster, RandomShp, sp=TRUE)
head(RandomDat.spdf1)
#Extract NLCD Data
RandomDat.spdf2 <- raster::extract(NLCDRaster, RandomDat.spdf1 , sp=TRUE)
head(RandomDat.spdf2)
# Convert back to data frame
RandomDat.df1 <- as.data.frame(RandomDat.spdf2, stringsAsFactors=FALSE)
head(RandomDat.df1)
# ReOrder last two columns
RandomDat.df <- RandomDat.df1[,c(1,6:7,2:5)]
head(RandomDat.df)
# Rename some columns
colnames(RandomDat.df)[1] <- c("Species")
colnames(RandomDat.df)[6:7] <- c("PopDen", "ClassNumber")
# Rename NEAR_DIST as RoadDist
names(RandomDat.df)[names(RandomDat.df) == "NEAR_DIST"] <- "RoadDist"
#
RandomDat.df$Species <- paste0(SpecShort, "Random")
nrow(RandomDat.df)
# Remove any no data
RandomDat.df <- na.omit(RandomDat.df)
# Drop unneeded columns
RandomDat.df <- RandomDat.df[,!(names(RandomDat.df) %in% c("NEAR_FID"))]
head(RandomDat.df)
## Merge ClassNumber values with Class Names into data frame
RandomDat.df2 <- merge(RandomDat.df, NLCDClassNames.df, by="ClassNumber")
head(RandomDat.df2)
# Delete class of "Open Water"
RandomDat.df2 <- RandomDat.df2[!RandomDat.df2$ClassName=="Open Water",]
# Save Random Data
write.table(RandomDat.df2, sep = ",", col.names=NA, file=paste0(Species, "_RandomPoints_Thinned1km_10000RoadProximityPopDensity.csv"))
###
#
####################################################
#
# Convert PresenceDat.df to spatial points data frame
head(PresenceDat.df)
PresenceDat.df[1:10,]
xy <- PresenceDat.df[,c("longitd","latitud")]
PresenceDat.spdf <- SpatialPointsDataFrame(coords = xy, data = PresenceDat.df, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
#
#str(PresenceDat.spdf)
#
PresenceDat.spdf@data[1:10,]
PresenceDat.spdf2 <- raster::extract(NLCDRaster, PresenceDat.spdf, sp=TRUE)
PresenceDat.spdf2@data[1:10,]
# Convert back to data frame
PresenceDatComb.df <- as.data.frame(PresenceDat.spdf2, stringsAsFactors=FALSE)
head(PresenceDatComb.df)
# Delete last two columns
PresenceDatComb.df <- PresenceDatComb.df[,-c((ncol(PresenceDatComb.df)-1):(ncol(PresenceDatComb.df)))]
# Rename last column to "ClassNumber"
colnames(PresenceDatComb.df)[ncol(PresenceDatComb.df)] <- "ClassNumber"
nrow(PresenceDatComb.df)
# Omit any rows with NA values in ClassNumber
PresenceDatComb.df <- PresenceDatComb.df[complete.cases(PresenceDatComb.df[,c("ClassNumber")]),]
nrow(PresenceDatComb.df)
#str(PresenceDat.df2)
################
## Merge ClassNumber values with Class Names into data frame
PresenceDatComb.df2 <- merge(PresenceDatComb.df, NLCDClassNames.df, by="ClassNumber")
head(PresenceDatComb.df2)
#str(PresenceDatComb.df2)
# Delete class of "Open Water"
nrow(PresenceDatComb.df2)
PresenceDatComb.df2 <- PresenceDatComb.df2[!PresenceDatComb.df2$ClassName=="Open Water",]
#
PresenceDatComb.df3 <- arrange(PresenceDatComb.df2, ClassName)
head(PresenceDatComb.df3)
tail(PresenceDatComb.df3)
nrow(PresenceDatComb.df3)
# Save data with original land cover classes
write.table(PresenceDatComb.df3, sep = ",", col.names=NA, file=paste0(Species, "_Thinned", SpatFiltBuff, Units,"_OriginalpresNLCD_", "1kPopDenRoaddat.csv"))
#
################################
######
## Read back in clipped and spatially thinned target group points
TargetGroupThinShp <- readOGR(".", paste0(Species, "_TargetGroupHiRes_Thinned", SpatFiltBuff,"kmrd"))  # NOTE: readShapePoly was giving error, so used readOGR instead
head(TargetGroupThinShp)
nrow(TargetGroupThinShp)
#
TargetGroupDat.spdf1 <- raster::extract(PopDenRaster, TargetGroupThinShp, sp=TRUE)
head(TargetGroupDat.spdf1)
TargetGroupDat.spdf2 <- raster::extract(NLCDRaster, TargetGroupDat.spdf1, sp=TRUE)
head(TargetGroupDat.spdf2)
# Convert back to data frame
TargetGroupDat.df <- as.data.frame(TargetGroupDat.spdf2, stringsAsFactors=FALSE)
head(TargetGroupDat.df)
# Delete last two columns
TargetGroupDat.df <- TargetGroupDat.df[,-c((ncol(TargetGroupDat.df)-1):(ncol(TargetGroupDat.df)))]
#
# Rename "sc2010pop1k" to "PopDen"
colnames(TargetGroupDat.df) <- gsub(PopGridName, "PopDen", colnames(TargetGroupDat.df))
colnames(TargetGroupDat.df) <- gsub(NLCDRasterName, "ClassNumber", colnames(TargetGroupDat.df))
names(TargetGroupDat.df)[names(TargetGroupDat.df) == "NEAR_DIST"] <- "RoadDist"
nrow(TargetGroupDat.df)
# Omit any rows with NA values for PopDen
# Make function to delete rows with NAs in a specified column
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}
TargetGroupDat.df <- completeFun(TargetGroupDat.df, "PopDen")
nrow(TargetGroupDat.df)
head(TargetGroupDat.df)
#
TargetGroupDat.df$RawCount <- 1
head(TargetGroupDat.df)
nrow(TargetGroupDat.df)
#
# Omit any rows with NA values in ClassNumber
TargetGroupDat.df <- TargetGroupDat.df[complete.cases(TargetGroupDat.df[,c("ClassNumber")]),]
nrow(TargetGroupDat.df)
#str(TargetGroupDat.df2)
################
## Merge ClassNumber values with Class Names into data frame
TargetGroupDat.df2 <- merge(TargetGroupDat.df, NLCDClassNames.df, by="ClassNumber")
head(TargetGroupDat.df2)
#str(TargetGroupDatComb.df2)
# Delete class of "Open Water"
nrow(TargetGroupDat.df2)
TargetGroupDat.df2 <- TargetGroupDat.df2[!TargetGroupDat.df2$ClassName=="Open Water",]
#
TargetGroupDat.df3 <- arrange(TargetGroupDat.df2, ClassName)
head(TargetGroupDat.df3)
tail(TargetGroupDat.df3)
nrow(TargetGroupDat.df3)
# Delete "CID" column
TargetGroupDat.df3$CID <- NULL
#
# Save data with original land cover classes
write.table(TargetGroupDat.df3, sep = ",", col.names=NA, file=paste0(Species, "_TargetGroupThinned", SpatFiltBuff, Units,"_OriginalpresNLCD_", "1kPopDenRoaddat.csv"))
#
##
###################################################################################################
## Take above csv file from PresenceDatComb.df3 and add new column RevClassName with revised
## class based on examination of photos and satellite imagery using iNaturalist url
## Drop records with no iNaturalist images
####################################################################################################
## Read in revised land cover data from iNaturalist 
setwd(InDirect)
PresenceDatComb.df3 <- data.frame(read.csv(paste0(Species, "_Thinned", SpatFiltBuff, Units,"_OriginalpresNLCD_", "1kPopDenRoaddatRev.csv"), stringsAsFactors=FALSE))
head(PresenceDatComb.df3)
nrow(PresenceDatComb.df3)
#str(PresenceDatComb.df4)
# Rename NEAR_DIST as RoadDist
names(PresenceDatComb.df3)[names(PresenceDatComb.df3) == "NEAR_DIST"] <- "RoadDist"
###
## Read in random and species target group data 
RandomDat.df2 <- data.frame(read.csv(paste0(Species, "_RandomPoints_Thinned1km_10000RoadProximityPopDensity.csv"), stringsAsFactors=FALSE))
head(RandomDat.df2)
#
TargetGroupDat.df3 <- data.frame(read.csv(paste0(Species, "_TargetGroupThinned", SpatFiltBuff, Units,"_OriginalpresNLCD_", "1kPopDenRoaddat.csv"), stringsAsFactors=FALSE))
head(TargetGroupDat.df3)
#
#####################
### Check correlation between PopDen and RoadDist for Random Data
PresenceDatComb.df3a <- subset(PresenceDatComb.df3, select = c("PopDen", "RoadDist"))
head(PresenceDatComb.df3a)
#
CorrRes1 <- corr.test(PresenceDatComb.df3a, method="spearman") # Holm's correction default
#str(CorrRes)
CorrRes1$r
CorrRes1$p
# Save above correlation output
out <- paste0("   ")
cat(out, file=paste0(Species, "PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out1 <-capture.output(print(CorrRes1$r))
cat(out1, file=paste0(Species, "PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out2 <- paste0("\n", "p-value with Holm Correction for Spearman Rank Order Correlations: ")
cat(out2, file=paste0(Species, "PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out3 <-capture.output(print(CorrRes1$p))
cat(out3, file=paste0(Species, "PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
##
#### Plot Data
p <- ggscatter(PresenceDatComb.df3a, x = "PopDen", y = "RoadDist",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.coef.size = 4, cor.method = "spearman",
          xlab = "Population Density per km", ylab = "Distance to Road (m)",
          title = paste0("                          ", LongSpecies))
#
p1 <- p +
  font("xlab", size = 18, face = "bold")+
  font("ylab", size = 18, face = "bold")+
  font("xy.text", size = 18, face = "bold")+
  font("title", face="bold")
#
#str(p1)
#
p2 <- p1 + theme(axis.title.x = element_text(vjust = -1)) # increased vertical space for xlab
p3 <- p2 + theme(axis.title.y = element_text(vjust = +1)) # increased vertical space for xlab
p4 <- p3 + theme(axis.line = element_line(colour = 'black', size = 1)) # increased line thickness
p5 <- p4 + theme(axis.ticks = element_line(colour = 'black', size = 1)) # increased tick thickness
#
ggarrange(p5, ncol=1) %>%
  ggexport(filename = paste0(InDirect, Species, "PopDenVsRoadDist.tiff"))
##
###############
### Check correlation between PopDen and RoadDist for Random Data
RandomDat.df2a <- subset(RandomDat.df2, select = c("PopDen", "RoadDist"))
head(RandomDat.df2a)
#
CorrRes1 <- corr.test(RandomDat.df2a, method="spearman") # Holm's correction default
#str(CorrRes)
CorrRes1$r
CorrRes1$p
# Save above correlation output
out <- paste0("   ")
cat(out, file=paste0(Species, "RandomDat_PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out1 <-capture.output(print(CorrRes1$r))
cat(out1, file=paste0(Species, "RandomDat_PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out2 <- paste0("\n", "p-value with Holm Correction for Spearman Rank Order Correlations: ")
cat(out2, file=paste0(Species, "RandomDat_PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out3 <-capture.output(print(CorrRes1$p))
cat(out3, file=paste0(Species, "RandomDat_PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
##
#### Plot Data
p <- ggscatter(RandomDat.df2a, x = "PopDen", y = "RoadDist",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.coef.size = 4, cor.method = "spearman",
          xlab = "Population Density per km", ylab = "Distance to Road (m)",
          title = paste0("            Random Data for ", LongSpecies))
#
p1 <- p +
  font("xlab", size = 18, face = "bold")+
  font("ylab", size = 18, face = "bold")+
  font("xy.text", size = 18, face = "bold")+
  font("title", face="bold")
#
#str(p1)
#
p2 <- p1 + theme(axis.title.x = element_text(vjust = -1)) # increased vertical space for xlab
p3 <- p2 + theme(axis.title.y = element_text(vjust = +1)) # increased vertical space for xlab
p4 <- p3 + theme(axis.line = element_line(colour = 'black', size = 1)) # increased line thickness
p5 <- p4 + theme(axis.ticks = element_line(colour = 'black', size = 1)) # increased tick thickness
#
ggarrange(p5, ncol=1) %>%
  ggexport(filename = paste0(InDirect, Species, "RandomDat_PopDenVsRoadDist.tiff"))
##
###############
### Check correlation between PopDen and RoadDist for Target Group Data
TargetGroupDat.df3a <- subset(TargetGroupDat.df3, select = c("PopDen", "RoadDist"))
head(TargetGroupDat.df3a)
#
CorrRes1 <- corr.test(TargetGroupDat.df3a, method="spearman") # Holm's correction default
#str(CorrRes)
CorrRes1$r
CorrRes1$p
# Save above correlation output
out <- paste0("   ")
cat(out, file=paste0(Species, "TargetGroup_PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out1 <-capture.output(print(CorrRes1$r))
cat(out1, file=paste0(Species, "TargetGroup_PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out2 <- paste0("\n", "p-value with Holm Correction for Spearman Rank Order Correlations: ")
cat(out2, file=paste0(Species, "TargetGroup_PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
out3 <-capture.output(print(CorrRes1$p))
cat(out3, file=paste0(Species, "TargetGroup_PopDenVsRoadDist_Correlation.txt"), sep="\n", append=TRUE)
##
#### Plot Data
p <- ggscatter(TargetGroupDat.df3a, x = "PopDen", y = "RoadDist",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.coef.size = 4, cor.method = "spearman",
          xlab = "Population Density per km", ylab = "Distance to Road (m)",
          title = paste0("        Target Group Data for ", LongSpecies))
#
p1 <- p +
  font("xlab", size = 18, face = "bold")+
  font("ylab", size = 18, face = "bold")+
  font("xy.text", size = 18, face = "bold")+
  font("title", face="bold")
#
#str(p1)
#
p2 <- p1 + theme(axis.title.x = element_text(vjust = -1)) # increased vertical space for xlab
p3 <- p2 + theme(axis.title.y = element_text(vjust = +1)) # increased vertical space for xlab
p4 <- p3 + theme(axis.line = element_line(colour = 'black', size = 1)) # increased line thickness
p5 <- p4 + theme(axis.ticks = element_line(colour = 'black', size = 1)) # increased tick thickness
#
ggarrange(p5, ncol=1) %>%
  ggexport(filename = paste0(InDirect, Species, "TargetGroupDat_PopDenVsRoadDist.tiff"))
##
##
##################################################################################################
### Histogram Bin Ratio Weighting (HBRWt) for correcting sample bias
### Calculate histograms of PopDen and RoadDist for random, target group, and species with same bins. Divide target group 
### percentage into random percentage for each bin to obtain a correction factor to weight
### the species count in the bin. Use zero weight, 50% weight and 100% weight to species and compare 
### results
################################################################################################### 
##  First normalize all values of PopDen across all three datasets
# Find the minimum and maximum values from all PopDen data to use 
# for normalization from zero to one
PopDenMin <- min(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, PresenceDatComb.df3$PopDen)
PopDenMax <- max(RandomDat.df2$PopDen, TargetGroupDat.df3$PopDen, PresenceDatComb.df3$PopDen)
# Calculate Normalized PopDen, PopDenN
RandomDat.df2$PopDenN <- (RandomDat.df2$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
head(RandomDat.df2)
TargetGroupDat.df3$PopDenN <- (TargetGroupDat.df3$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
PresenceDatComb.df3$PopDenN <- (PresenceDatComb.df3$PopDen - PopDenMin)/(PopDenMax - PopDenMin)
## Prepare histogram for random points PopDen
png(file = paste0(Species, "_RandomPoints_PopDenNHistogram.png"))
randomhistres <- hist(RandomDat.df2$PopDenN, breaks=seq(0,1,by=0.001))
plot(randomhistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
## Prepare histogram for target group PopDen
png(file = paste0(Species, "_TargetGroup_PopDenNHistogram.png"))
targethistres <- hist(TargetGroupDat.df3$PopDenN, breaks=seq(0,1,by=0.001))
plot(targethistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
## Prepare histogram for presence group PopDen
png(file = paste0(Species, "_Presence_PopDenNHistogram.png"))
presencehistres <- hist(PresenceDatComb.df3$PopDenN, breaks=seq(0,1,by=0.001))
plot(presencehistres, col="blue", border= "blue", xlab = "Normalized Population Density")
dev.off()
############
length(targethistres$count)
randomhistres$count
presencehistres$count

## Assemble a dataframe with all histogram results
##
PopDenHBRWt.df <- data.frame(cbind(randomhistres$count, targethistres$count, presencehistres$count))
head(PopDenHBRWt.df)
tail(PopDenHBRWt.df)
colnames(PopDenHBRWt.df) <- c("Random", "TargetGroup", "Presence")
# Calculate frequency percent for each column (column will sum to 1)
PopDenHBRWt.df$RandomFreqPerc <- PopDenHBRWt.df$Random/sum(PopDenHBRWt.df$Random) 
PopDenHBRWt.df$TargetGroupFreqPerc <- PopDenHBRWt.df$TargetGroup/sum(PopDenHBRWt.df$TargetGroup) 
PopDenHBRWt.df$PresenceFreqPerc <- PopDenHBRWt.df$Presence/sum(PopDenHBRWt.df$Presence)
# Replace all values of zero with 50% of minimum value above zero
PopDenHBRWt.df$RandomFreqPerc[PopDenHBRWt.df$RandomFreqPerc==0] <- unique(sort(PopDenHBRWt.df$RandomFreqPerc))[2]*0.5
PopDenHBRWt.df$TargetGroupFreqPerc[PopDenHBRWt.df$TargetGroupFreqPerc==0] <- unique(sort(PopDenHBRWt.df$TargetGroupFreqPerc))[2]*0.5
PopDenHBRWt.df$PresenceFreqPerc[PopDenHBRWt.df$PresenceFreqPerc==0] <- unique(sort(PopDenHBRWt.df$PresenceFreqPerc))[2]*0.5
# Calculate the PopDenHBRWt
PopDenHBRWt.df$PopDenHBRWt <- PopDenHBRWt.df$RandomFreqPerc/PopDenHBRWt.df$TargetGroupFreqPerc
# Replace all values of zero with 50% of minimum value above zero
PopDenHBRWt.df$PopDenHBRWt[PopDenHBRWt.df$PopDenHBRWt==0] <- unique(sort(PopDenHBRWt.df$PopDenHBRWt))[2]*0.5
# Provide bin values
PopDenHBRWt.df$Bin <- seq(0.001,1,by=0.001)
head(PopDenHBRWt.df)
## Assign Presence data PopDenN values to bins corresponding to PopDenHBRWt.df$Bin
PresenceDatComb.df3$PopDenNBinCode <- .bincode(PresenceDatComb.df3$PopDenN, seq(0,1,by=0.001), include.lowest=TRUE)
head(PresenceDatComb.df3)
# Retrieve HBRWt for PopDenN
PresenceDatComb.df3$PopDenNHBRWt <- PopDenHBRWt.df[PresenceDatComb.df3$PopDenNBinCode,7]
head(PresenceDatComb.df3)
PopDenHBRWt.df[1:20,]
####################################
### Calculate weighting for RoadDist using new PopDenN bins from random data
##
## First normalize all values of RoadDist across all three datasets
# Find the minimum and maximum values from all PopDen data to use 
# for normalization from zero to one
RoadDistMin <- min(RandomDat.df2$RoadDist, TargetGroupDat.df3$RoadDist, PresenceDatComb.df3$RoadDist)
RoadDistMax <- max(RandomDat.df2$RoadDist, TargetGroupDat.df3$RoadDist, PresenceDatComb.df3$RoadDist)
# Calculate Normalized RoadDist, RoadDistN
RandomDat.df2$RoadDistN <- (RandomDat.df2$RoadDist - RoadDistMin)/(RoadDistMax - RoadDistMin)
head(RandomDat.df2)
TargetGroupDat.df3$RoadDistN <- (TargetGroupDat.df3$RoadDist - RoadDistMin)/(RoadDistMax - RoadDistMin)
PresenceDatComb.df3$RoadDistN <- (PresenceDatComb.df3$RoadDist - RoadDistMin)/(RoadDistMax - RoadDistMin)
###
# Find breaks for values of PopDenN that will be considered for nested histograms of RoadDistN
# First calculate cumulative frequencies of PopDenN
PopDenHBRWt.df$RandomCumSum <- cumsum(PopDenHBRWt.df$Random)
RandomCumSumDiff <- diff(PopDenHBRWt.df$RandomCumSum)
PopDenHBRWt.df$RandomCumSumDiff <- c(PopDenHBRWt.df$Random[1], RandomCumSumDiff)
#
sum(PopDenHBRWt.df$Random)
sum(PopDenHBRWt.df$TargetGroup)
sum(PopDenHBRWt.df$Presence)
max(PopDenHBRWt.df$Bin)
## Assign bin values for every PopDenN bin with at least 100 members
PopDenHBRWt.df$PopDenNBin100 <- 999
nrow(PopDenHBRWt.df)
#for(i in 1:nrow(PopDenHBRWt.df)){
for(i in 1:20) {
  #i=5
  if(PopDenHBRWt.df[i,(ncol(PopDenHBRWt.df)-1)]>99) {
    PopDenHBRWt.df[i,ncol(PopDenHBRWt.df)] <- PopDenHBRWt.df[i,(ncol(PopDenHBRWt.df)-3)]
  } else {
    sumRandomCumSumDiff <- 0
    for(j in 1:i) {
      #j=1
      if(PopDenHBRWt.df[i-j,(ncol(PopDenHBRWt.df))]==999) {
        sumRandomCumSumDiff <- sumRandomCumSumDiff + PopDenHBRWt.df[i-j,(ncol(PopDenHBRWt.df)-1)]
        if(PopDenHBRWt.df[i,(ncol(PopDenHBRWt.df)-1)] + sumRandomCumSumDiff > 99) {
          PopDenHBRWt.df[i,ncol(PopDenHBRWt.df)] <- PopDenHBRWt.df[i,(ncol(PopDenHBRWt.df)-3)]
        } 
      } else {
        break
      }
    }
  }
}
#
PopDenNRoadDistNBreaks <- PopDenHBRWt.df$PopDenNBin100
PopDenNRoadDistNBreaks <- PopDenNRoadDistNBreaks[PopDenNRoadDistNBreaks<999]
PopDenNRoadDistNBreaks <- c(0,PopDenNRoadDistNBreaks, 1)
#
# Group all three datasets according to above bin breaks
# Then perform RoadDistN histogram on each group and calculate weighting value within groups
RandomDatGroups.list <- split(RandomDat.df2, cut(RandomDat.df2$PopDenN, PopDenNRoadDistNBreaks, include.lowest=TRUE))
TargetGroupDatGroups.list <- split(TargetGroupDat.df3, cut(TargetGroupDat.df3$PopDenN, PopDenNRoadDistNBreaks, include.lowest=TRUE))
PresenceDatGroups.list <- split(PresenceDatComb.df3, cut(PresenceDatComb.df3$PopDenN, PopDenNRoadDistNBreaks, include.lowest=TRUE))
#str(RandomDatGroups.list)
nrow(RandomDat.df2)
nrow(TargetGroupDat.df3)
nrow(PresenceDatComb.df3)
#
Randomlengthlist <- list()
TargetGrouplengthlist <- list()
Presencelengthlist <- list()
names(RandomDatGroups.list)[[1]]
#RandomDatGroups.list[[k]]$PopDenN
for(k in 1:length(RandomDatGroups.list)) {
  #k=1
  Randomlengthlist[[k]] <- nrow(RandomDatGroups.list[[k]])
  TargetGrouplengthlist[[k]] <- nrow(TargetGroupDatGroups.list[[k]])
  Presencelengthlist[[k]] <- nrow(PresenceDatGroups.list[[k]])
  ## Prepare histogram for random points PopDen
  png(file = paste0(Species, "_RandomPoints_", names(RandomDatGroups.list)[[k]], "_RoadDistNHistogram.png"))
  randomhistres <- hist(RandomDatGroups.list[[k]]$RoadDistN, breaks=seq(0,1,by=0.001))
  plot(randomhistres, col="blue", border= "blue", xlab = "Normalized Distance to Road")
  dev.off()
  ## Prepare histogram for target group PopDen
  png(file = paste0(Species, "_TargetGroup_", names(TargetGroupDatGroups.list)[[k]], "_RoadDistNHistogram.png"))
  targethistres <- hist(TargetGroupDatGroups.list[[k]]$RoadDistN, breaks=seq(0,1,by=0.001))
  plot(targethistres, col="blue", border= "blue", xlab = "Normalized Distance to Road")
  dev.off()
  ## Prepare histogram for Presence Group
  png(file = paste0(Species, "_Presence_", names(PresenceDatGroups.list)[[k]], "_RoadDistNHistogram.png"))
  presencehistres <- hist(PresenceDatGroups.list[[k]]$RoadDistN, breaks=seq(0,1,by=0.001))
  plot(presencehistres, col="blue", border= "blue", xlab = "Normalized Distance to Road")
  dev.off()
  ############
  length(targethistres$count)
  randomhistres$count
  presencehistres$count
  #
  ## Assemble a dataframe with all histogram results
  ##
  RoadDistHBRWt.df <- data.frame(cbind(randomhistres$count, targethistres$count, presencehistres$count))
  head(RoadDistHBRWt.df)
  tail(RoadDistHBRWt.df)
  colnames(RoadDistHBRWt.df) <- c("Random", "TargetGroup", "Presence")
  # Calculate frequency percent for each column (column will sum to 1)
  RoadDistHBRWt.df$RandomFreqPerc <- RoadDistHBRWt.df$Random/sum(RoadDistHBRWt.df$Random) 
  RoadDistHBRWt.df$TargetGroupFreqPerc <- RoadDistHBRWt.df$TargetGroup/sum(RoadDistHBRWt.df$TargetGroup) 
  RoadDistHBRWt.df$PresenceFreqPerc <- RoadDistHBRWt.df$Presence/sum(RoadDistHBRWt.df$Presence)
  # Replace all values of zero with 50% of minimum value above zero
  RoadDistHBRWt.df$RandomFreqPerc[RoadDistHBRWt.df$RandomFreqPerc==0] <- unique(sort(RoadDistHBRWt.df$RandomFreqPerc))[2]*0.5
  RoadDistHBRWt.df$TargetGroupFreqPerc[RoadDistHBRWt.df$TargetGroupFreqPerc==0] <- unique(sort(RoadDistHBRWt.df$TargetGroupFreqPerc))[2]*0.5
  RoadDistHBRWt.df$PresenceFreqPerc[RoadDistHBRWt.df$PresenceFreqPerc==0] <- unique(sort(RoadDistHBRWt.df$PresenceFreqPerc))[2]*0.5
  # Calculate the RoadDistHBRWt
  RoadDistHBRWt.df$RoadDistHBRWt <- RoadDistHBRWt.df$RandomFreqPerc/RoadDistHBRWt.df$TargetGroupFreqPerc
  # Replace all values of zero with 50% of minimum value above zero
  RoadDistHBRWt.df$RoadDistHBRWt[RoadDistHBRWt.df$RoadDistHBRWt==0] <- unique(sort(RoadDistHBRWt.df$RoadDistHBRWt))[2]*0.5
  # Provide bin values
  RoadDistHBRWt.df$Bin <- seq(0.001,1,by=0.001)
  head(RoadDistHBRWt.df)
  ## Assign Presence data RoadDistN values to bins corresponding to RoadDistHBRWt.df$Bin
  PresenceDatGroups.list[[k]]$RoadDistBinCode <- .bincode(PresenceDatGroups.list[[k]]$RoadDistN, seq(0,1,by=0.001), include.lowest=TRUE)
  head(PresenceDatGroups.list[[k]])
  # Retrieve HBRWt for RoadDistN
  PresenceDatGroups.list[[k]]$RoadDistNHBRWt <- RoadDistHBRWt.df[PresenceDatGroups.list[[k]]$RoadDistBinCode,7]
  head(PresenceDatGroups.list[[k]])
}
########
sum(unlist(Presencelengthlist))
sum(unlist(TargetGrouplengthlist))
sum(unlist(Randomlengthlist))
#
PresenceDatComb.df4 <- do.call(rbind, PresenceDatGroups.list)
head(PresenceDatComb.df4)
tail(PresenceDatComb.df4)
nrow(PresenceDatComb.df4)
############################################################################################
## Calculate combination PopDenNHBRWt and RoadDistNHBRWt 
PresenceDatComb.df4$PopRoadHBRWt <- PresenceDatComb.df4$PopDenNHBRWt * PresenceDatComb.df4$RoadDistNHBRWt
# Calculate half distance weights for all HBRWts
PresenceDatComb.df4$PopDenNHBRWtH <- NA
for(i in 1:(nrow(PresenceDatComb.df4))) {
  #i=1
  if(PresenceDatComb.df4$PopDenNHBRWt[i] < 1) {
    PresenceDatComb.df4$PopDenNHBRWtH[i] <- (0.5*(1 - PresenceDatComb.df4$PopDenNHBRWt[i]))+ PresenceDatComb.df4$PopDenNHBRWt[i]
  } else if (PresenceDatComb.df4$PopDenNHBRWt[i] > 1) {
    PresenceDatComb.df4$PopDenNHBRWtH[i] <- PresenceDatComb.df4$PopDenNHBRWt[i] - (0.5*(PresenceDatComb.df4$PopDenNHBRWt[i] - 1)) 
  } else {
    PresenceDatComb.df4$PopDenNHBRWtH[i] <- PresenceDatComb.df4$PopDenNHBRWt[i]
  }
}
##
PresenceDatComb.df4$RoadDistNHBRWtH <- NA
for(i in 1:(nrow(PresenceDatComb.df4))) {
  #i=1
  if(PresenceDatComb.df4$RoadDistNHBRWt[i] < 1) {
    PresenceDatComb.df4$RoadDistNHBRWtH[i] <- (0.5*(1 - PresenceDatComb.df4$RoadDistNHBRWt[i]))+ PresenceDatComb.df4$RoadDistNHBRWt[i]
  } else if (PresenceDatComb.df4$RoadDistNHBRWt[i] > 1) {
    PresenceDatComb.df4$RoadDistNHBRWtH[i] <- PresenceDatComb.df4$RoadDistNHBRWt[i] - (0.5*(PresenceDatComb.df4$RoadDistNHBRWt[i] - 1)) 
  } else {
    PresenceDatComb.df4$RoadDistNHBRWtH[i] <- PresenceDatComb.df4$RoadDistNHBRWt[i]
  }
}
##
PresenceDatComb.df4$PopRoadHBRWtH <- NA
for(i in 1:(nrow(PresenceDatComb.df4))) {
  #i=1
  if(PresenceDatComb.df4$PopRoadHBRWt[i] < 1) {
    PresenceDatComb.df4$PopRoadHBRWtH[i] <- (0.5*(1 - PresenceDatComb.df4$PopRoadHBRWt[i]))+ PresenceDatComb.df4$PopRoadHBRWt[i]
  } else if (PresenceDatComb.df4$PopRoadHBRWt[i] > 1) {
    PresenceDatComb.df4$PopRoadHBRWtH[i] <- PresenceDatComb.df4$PopRoadHBRWt[i] - (0.5*(PresenceDatComb.df4$PopRoadHBRWt[i] - 1)) 
  } else {
    PresenceDatComb.df4$PopRoadHBRWtH[i] <- PresenceDatComb.df4$PopRoadHBRWt[i]
  }
}
#################################
# Normalize all weighting factors to sum to total raw count for calculating expected counts
TotRawCount <- sum(PresenceDatComb.df4$RawCount)
TotPopDenNHBRWt <- sum(PresenceDatComb.df4$PopDenNHBRWt)
TotPopDenNHBRWtH <- sum(PresenceDatComb.df4$PopDenNHBRWtH)
TotRoadDistNHBRWt <- sum(PresenceDatComb.df4$RoadDistNHBRWt)
TotRoadDistNHBRWtH <- sum(PresenceDatComb.df4$RoadDistNHBRWtH)
TotPopRoadHBRWt <- sum(PresenceDatComb.df4$PopRoadHBRWt)
TotPopRoadHBRWtH <- sum(PresenceDatComb.df4$PopRoadHBRWtH)
# Normalize weights to sum to TotRawCount
PresenceDatComb.df4$PopDenNHBRWtN <- PresenceDatComb.df4$PopDenNHBRWt/(TotPopDenNHBRWt/TotRawCount)
PresenceDatComb.df4$PopDenNHBRWtHN <- PresenceDatComb.df4$PopDenNHBRWtH/(TotPopDenNHBRWtH/TotRawCount)
PresenceDatComb.df4$RoadDistNHBRWtN <- PresenceDatComb.df4$RoadDistNHBRWt/(TotRoadDistNHBRWt/TotRawCount)
PresenceDatComb.df4$RoadDistNHBRWtHN <- PresenceDatComb.df4$RoadDistNHBRWtH/(TotRoadDistNHBRWtH/TotRawCount)
PresenceDatComb.df4$PopRoadHBRWtN <- PresenceDatComb.df4$PopRoadHBRWt/(TotPopRoadHBRWt/TotRawCount)
PresenceDatComb.df4$PopRoadHBRWtHN <- PresenceDatComb.df4$PopRoadHBRWtH/(TotPopRoadHBRWtH/TotRawCount)
nrow(PresenceDatComb.df4)
#
TotPopDenNHBRWtN <- sum(PresenceDatComb.df4$PopDenNHBRWtN)
TotPopDenNHBRWtHN <- sum(PresenceDatComb.df4$PopDenNHBRWtHN)
TotRoadDistNHBRWtN <- sum(PresenceDatComb.df4$RoadDistNHBRWtN)
TotRoadDistNHBRWtHN <- sum(PresenceDatComb.df4$RoadDistNHBRWtHN)
TotPopRoadHBRWtN <- sum(PresenceDatComb.df4$PopRoadHBRWtN)
TotPopRoadHBRWtHN <- sum(PresenceDatComb.df4$PopRoadHBRWtHN)
#
##########################################
## Read in tabulated random availability of land cover classes 
## obtained from NLCD raster clipped to MaxEnt core convex hull polygon shapefile
##
RandomNLCD.df <- as.data.frame(read.dbf(NLCDMxtCoreTableName), stringsAsFactors=FALSE)
#str(RandomNLCD.df)
# Delete first column
RandomNLCD.df <- RandomNLCD.df[,c(-1)]
# Fix column names
colClean <- function(x){ colnames(x) <- gsub("VALUE_", "", colnames(x)); x }
RandomNLCD.df <- colClean(RandomNLCD.df)
# Transpose data frame
RandomNLCD.df1 <- as.data.frame(t(RandomNLCD.df), stringsAsFactors=FALSE)
#str(RandomNLCD.df1)
# Use row names as column "ClassNumber"
setDT(RandomNLCD.df1, keep.rownames = "ClassNumber")
colnames(RandomNLCD.df1)[2] <- "Count"
RandomNLCD.df1$ClassNumber <- as.integer(RandomNLCD.df1$ClassNumber)
# Re-associate with land cover class names
## Merge ClassNumber values with Class Names into data frame
RandomNLCD.df2 <- merge(RandomNLCD.df1, NLCDClassNames.df, by.x="ClassNumber", by.y="ClassNumber")
head(RandomNLCD.df2)
#str(RandomDat.df3)
# Delete open water
RandomNLCD.df2 <- RandomNLCD.df2[!RandomNLCD.df2$ClassName=="Open Water",]
#str(RandomNLCD.df2)
head(RandomNLCD.df2)
# Calculate RandomPercAvail
RandomNLCD.df2$RandomPercAvail <- RandomNLCD.df2$Count/sum(RandomNLCD.df2$Count)
# Calculate expected count, which is RandomPercAvail * Total Actual PopDenInd across all classes except water
RandomNLCD.df2$ExpectedCount <-  RandomNLCD.df2$RandomPercAvail*TotRawCount
TotExpectedCount <- sum(RandomNLCD.df2$ExpectedCount)
## Identify and combine classes that have expected counts below 5 (c.f. Bryers et al. 1984, Manly et al. 2002)
CombineDropClasses.df <- RandomNLCD.df2[ which(RandomNLCD.df2$ExpectedCount<5),]
print(CombineDropClasses.df)
#str(DropClasses.df)
# Save Dropped Classes
write.table(CombineDropClasses.df, sep= ",", col.names=NA, file=paste0(Species, "CombinedNLCDLandCoverClassesforLandCoverUtilIndex.csv"))
#
# Save above list of dropped classes to text file
out1 <- paste0("List of Class Names combined (Open Water dropped) that had expected count less than 5")
cat(out1, file=paste0(Species, "ListCombinedNLCDLandCoverClassesforLandCoverUtiIndex.txt"), sep="\n", append=TRUE)
out2 <- capture.output(print(CombineDropClasses.df))
cat(out2, file=paste0(Species, "ListCombinedNLCDLandCoverClassesforLandCoverUtiIndex.txt"), sep="\n", append=TRUE)
#
CombineDropClasses <- CombineDropClasses.df$ClassName
## Examine data to determine which similar classes to combine and change names for combined classes
head(PresenceDatComb.df4)
PresenceDatComb.df4$NewClassName <- PresenceDatComb.df4$RevClassName
PresenceDatComb.df4a <- PresenceDatComb.df4
#PresenceDatComb.df4 <- PresenceDatComb.df4a
PresenceDatComb.df4$NewClassName
## Combine Mixed Forest and Evergreen Forest
#PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
#  {gsub("Mixed Forest", "Moxed Forest", x)}))
#  PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
#  {gsub("Evergreen Forest", "Evergreen/Mixed Forest", x)}))
#PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
#  {gsub("Moxed Forest", "Evergreen/Mixed Forest", x)}))
# Combine Developed, Medium Intensity with Developed, High Intensity
PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
  {gsub("Developed, Medium Intensity", "Developed, Med/Hi Intensity", x)}))
PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
  {gsub("Developed, High Intensity", "Developed, Med/Hi Intensity", x)}))
# Combined Barren, Woody Wetlands, Emergent Herbaceous, and Mixed Forest into Other
PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
  {gsub("Barren Land", "Other", x)}))
PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
  {gsub("Woody Wetlands", "Other", x)}))
PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
  {gsub("Emergent Herbaceous Wetlands", "Other", x)}))
PresenceDatComb.df4$NewClassName <- unlist(lapply(PresenceDatComb.df4$NewClassName, function(x)
  {gsub("Mixed Forest", "Other", x)}))
#
unique(PresenceDatComb.df4$NewClassName)
#
# Save data
write.table(PresenceDatComb.df4, sep = ",", col.names=NA, file=paste0(Species, "_Thinned", SpatFiltBuff, Units,"_presNLCD_", "1kPopDenIndexdat.csv"))
#
#### Combine various classes for RandomNLCD.df2
RandomNLCD.df2$NewClassName <- RandomNLCD.df2$ClassName
unique(RandomNLCD.df2$NewClassName)
## Combine Mixed Forest and Evergreen Forest
#RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
#  {gsub("Mixed Forest", "Moxed Forest", x)}))
#  RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
#  {gsub("Evergreen Forest", "Evergreen/Mixed Forest", x)}))
#RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
#  {gsub("Moxed Forest", "Evergreen/Mixed Forest", x)}))
# Combine Developed, Medium Intensity with Developed, High Intensity
RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
  {gsub("Developed, Medium Intensity", "Developed, Med/Hi Intensity", x)}))
RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
  {gsub("Developed, High Intensity", "Developed, Med/Hi Intensity", x)}))
# Combined Barren, Woody Wetlands, and Emergent Herbaceous into Other
RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
  {gsub("Barren Land", "Other", x)}))
RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
  {gsub("Woody Wetlands", "Other", x)}))
RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
  {gsub("Emergent Herbaceous Wetlands", "Other", x)}))
RandomNLCD.df2$NewClassName <- unlist(lapply(RandomNLCD.df2$NewClassName, function(x)
  {gsub("Mixed Forest", "Other", x)}))
#
unique(RandomNLCD.df2$NewClassName)
#  Aggregate count by NewClassName and recalculate RandomPercAvail and ExpectedCount
RandomNLCD.df5 <- aggregate(Count ~ NewClassName, data=RandomNLCD.df2, FUN=function(x) c(CountSum = sum(x)))
# Calculate RandomPercAvail
RandomNLCD.df5$RandomPercAvail <- RandomNLCD.df5$Count/sum(RandomNLCD.df5$Count)
# Calculate expected count, which is RandomPercAvail * Total Actual PopDenInd across all classes except water
RandomNLCD.df5$ExpectedCount <-  RandomNLCD.df5$RandomPercAvail*TotRawCount
TotExpectedCount2 <- sum(RandomNLCD.df5$ExpectedCount)
TotRandomPercAvail <- sum(RandomNLCD.df5$RandomPercAvail)
##
# Save set of all random data with dropped classes
## Save data
write.table(RandomNLCD.df5, sep = ",", col.names=NA, file=paste0(Species, "_HiRes_NLCD", Year, "LandCoverStatistics.csv"))
#####################################################################################################


###################
#
### Save names of classes in RandomNLCD.df5 not present in PresenceDatComb.df3,
## First collect classes in RandomNLCD.df5
RandomClasses <- RandomNLCD.df5$NewClassName
## Next, collect all classes in PresenceDatComb.df3
PresenceClasses <- unique(PresenceDatComb.df4$NewClassName)
## Find classes in random but not presence
NewClasses.df <- as.data.frame(setdiff(RandomClasses, PresenceClasses), stringsAsFactors=FALSE)
colnames(NewClasses.df) <- "ClassName"
# Save New Classes
write.table(NewClasses.df, sep= ",", col.names=NA, file=paste0(Species, "AddedRandomNLCDLandCoverClasses.csv"))
#
############################
## Obtain counts by Class Name
ActualDatTab.dfa <- aggregate(RawCount ~ NewClassName, data=PresenceDatComb.df4, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dfb <- aggregate(PopDenNHBRWtN ~ NewClassName, data=PresenceDatComb.df4, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dfc <- aggregate(PopDenNHBRWtHN ~ NewClassName, data=PresenceDatComb.df4, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dfd <- aggregate(RoadDistNHBRWtN ~ NewClassName, data=PresenceDatComb.df4, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dfe <- aggregate(RoadDistNHBRWtHN ~ NewClassName, data=PresenceDatComb.df4, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dff <- aggregate(PopRoadHBRWtN ~ NewClassName, data=PresenceDatComb.df4, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dfg <- aggregate(PopRoadHBRWtHN ~ NewClassName, data=PresenceDatComb.df4, FUN=function(x) c(Sum = sum(x)))
ActualDatTab.dfh <- Reduce(function(x,y) merge(x,y, by="NewClassName", all=TRUE), list(ActualDatTab.dfa, ActualDatTab.dfb, ActualDatTab.dfc, ActualDatTab.dfd, ActualDatTab.dfe, ActualDatTab.dff, ActualDatTab.dfg))
#
sum(ActualDatTab.dfh$RawCount)
sum(ActualDatTab.dfh$PopRoadHBRWtHN)
# Save raw data for original classes
write.table(ActualDatTab.dfh, sep = ",", col.names=NA, file=paste0(Species, "_Thinned", SpatFiltBuff, Units,"_LandClassDataCounts.csv"))
#
############################
# Read in data if needed
ActualDatTab.dfh <- data.frame(read.csv(paste0(Species, "_Thinned", SpatFiltBuff, Units,"_LandClassDataCounts.csv"), row.names=1, stringsAsFactors=FALSE))
head(ActualDatTab.dfh)
RandomNLCD.df5 <- data.frame(read.csv(paste0(Species, "_HiRes_NLCD", Year, "LandCoverStatistics.csv"), row.names=1, stringsAsFactors=FALSE))
head(RandomNLCD.df5)
#
ActualDatTab.dfi <- merge(ActualDatTab.dfh, RandomNLCD.df5, by="NewClassName", all=TRUE)
# Replace NA with zero
ActualDatTab.dfi[is.na(ActualDatTab.dfi)] <- 0
sum(ActualDatTab.dfi$RandomPercAvail)
ActualDatTab.dfi$RawPercUsed <- ActualDatTab.dfi$RawCount/(sum(ActualDatTab.dfi$RawCount))
ActualDatTab.dfi$PopDenNHBRWtNPercUsed <- ActualDatTab.dfi$PopDenNHBRWtN/(sum(ActualDatTab.dfi$RawCount))
ActualDatTab.dfi$PopDenNHBRWtHNPercUsed <- ActualDatTab.dfi$PopDenNHBRWtHN/(sum(ActualDatTab.dfi$RawCount))
ActualDatTab.dfi$RoadDistNHBRWtNPercUsed <- ActualDatTab.dfi$RoadDistNHBRWtN/(sum(ActualDatTab.dfi$RawCount))
ActualDatTab.dfi$RoadDistNHBRWtHNPercUsed <- ActualDatTab.dfi$RoadDistNHBRWtHN/(sum(ActualDatTab.dfi$RawCount))
ActualDatTab.dfi$PopRoadHBRWtNPercUsed <- ActualDatTab.dfi$PopRoadHBRWtN/(sum(ActualDatTab.dfi$RawCount))
ActualDatTab.dfi$PopRoadHBRWtHNPercUsed <- ActualDatTab.dfi$PopRoadHBRWtHN/(sum(ActualDatTab.dfi$RawCount))
## Save data
write.table(ActualDatTab.dfi, sep = ",", col.names=NA, file=paste0(Species, "_HiResNLCDPercentAvailableandUsed.csv"))
###############################################################################################################
## Calculate Bonferroni Confidence Intervals for whether RandomPercAvail falls within BCI for PercUsed
## based on raw data and various (six) bias adjustments for human populations and road proximity
## following Bryers et al. 1984 and Manly et al. 2002
###############################################################################################################
# Read back in data
ActualDatTab.dfi <- data.frame(read.csv(paste0(Species, "_HiResNLCDPercentAvailableandUsed.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
## Delete row names column
ActualDatTab.dfi <- ActualDatTab.dfi[,-1]
#
##
sum(ActualDatTab.dfi$ExpectedCount)
sum(ActualDatTab.dfi$PopRoadHBRWtHN)
sum(ActualDatTab.dfi$RandomPercAvail)
NumberComparisons <- nrow(ActualDatTab.dfi)
#
## Loop through the test for the raw and six bias adjustments for percent land cover used
## and calculate electivity indices
UsedColumnList <- seq(2, 8, 1)
UsedCategoryList <- colnames(ActualDatTab.dfi)[UsedColumnList]
PercUsedColumnList <- seq(12,18,1)
PercUsedCategoryList <- colnames(ActualDatTab.dfi)[PercUsedColumnList]
count=0
for(UsedCol in UsedColumnList) {
  #UsedCol=2
  PercUsedCol <- UsedCol + 10
  count=count+1
  ChiSqGFT_Results <- chisq.test(ActualDatTab.dfi[,UsedCol], p=ActualDatTab.dfi$RandomPercAvail)
  #str(ChiSqGFT_Results)
  #ChiSqGFT_Results$expected
  # Save above test results
  out1 <- paste0("Chi-Square Goodness of Fit Test for Differeces in Expected (Random) Versus ", UsedCategoryList[count]," Utilization of Land Classes")
  cat(out1, file=paste0(Species, "2016MilkweedLandCoverChiSquareGoodnessFitTestfor", UsedCategoryList[count],".txt"), sep="\n")
  out2 <- capture.output(print(ChiSqGFT_Results))
  cat(out2, file=paste0(Species, "2016MilkweedLandCoverChiSquareGoodnessFitTestfor", UsedCategoryList[count],".txt"), sep="\n", append=TRUE)
  #
  # Assuming above ChiSquare Test is significant, find bonferoni critical z-value of 0.05 with NumberComparisons categories for two-tailed 95% CI
  pc <- 0.05/(NumberComparisons)
  za <- qnorm(1-(pc/2))
  ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"BCI")] <- za*sqrt((ActualDatTab.dfi[,PercUsedCol]*(1-ActualDatTab.dfi[,PercUsedCol]))/sum(ActualDatTab.dfi$ExpectedCount))
  ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"LowerBCI")] <- ActualDatTab.dfi[,PercUsedCol] - ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"BCI")]
  ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"UpperBCI")] <- ActualDatTab.dfi[,PercUsedCol] + ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"BCI")]
  # If RandomPercAvail outside BCI, identify as Significant difference
  if(ChiSqGFT_Results$p.value<0.05) {
    ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"BCISig")] <- ifelse(ActualDatTab.dfi$RandomPercAvail <= ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"UpperBCI")] & ActualDatTab.dfi$RandomPercAvail >=  ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"LowerBCI")],  "NonSig",  "Significant")
    } else {
    ActualDatTab.dfi[, paste0(PercUsedCategoryList[count],"BCISig")] <- "NonSig"
    #
  }
  ActualDatTab.dfi[, paste0("ForateRat", PercUsedCategoryList[count])] <- ivlev_forage(ActualDatTab.dfi[,PercUsedCol], ActualDatTab.dfi$RandomPercAvail, log10=FALSE)
  ActualDatTab.dfi[, paste0("VS_Elect", PercUsedCategoryList[count])] <- vs_electivity(ActualDatTab.dfi[,PercUsedCol], ActualDatTab.dfi$RandomPercAvail)
  ActualDatTab.dfi[, paste0("VS_Elect", PercUsedCategoryList[count], "_10K")] <- Round2(ActualDatTab.dfi[, paste0("VS_Elect", PercUsedCategoryList[count])]*10000,0)
}
## Calculated Normalized VS_Elect
VS_ElectMin <- min(ActualDatTab.dfi$VS_ElectPopRoadHBRWtNPercUsed)
VS_ElectMax <- max(ActualDatTab.dfi$VS_ElectPopRoadHBRWtNPercUsed)
# Calculate Normalized PopDen, PopDenN
ActualDatTab.dfi$VS_ElectPopRoadHBRWtNPercUsedN_10K <- 10000*((ActualDatTab.dfi$VS_ElectPopRoadHBRWtNPercUsed - VS_ElectMin)/(VS_ElectMax - VS_ElectMin))
head(ActualDatTab.dfi)
# Calculate Percent of Total Normalized VS_Elect
ActualDatTab.dfi$VS_ElectPopRoadHBRWtNPercUsedNP <- ActualDatTab.dfi$VS_ElectPopRoadHBRWtNPercUsedN_10K/sum(ActualDatTab.dfi$VS_ElectPopRoadHBRWtNPercUsedN_10K)
#
## Save data
write.table(ActualDatTab.dfi, sep = ",", col.names=NA, file=paste0(Species, "_HiResNLCDActualPercentUsedStatsIndices.csv"))

########################################################################################
##  Summarize counts of description types per land classs
########################################################################################
## Read in High Resolution data with Road proximity corrected Population Index values, PopDenRoadInd
setwd(InDirect)
PresenceDatComb.df4 <- data.frame(read.csv(paste0(Species, "_Thinned", SpatFiltBuff, Units,"_presNLCD_", "1kPopDenIndexdat.csv"), stringsAsFactors=FALSE))
#
head(PresenceDatComb.df4)
sum(PresenceDatComb.df4$PopRoadHBRWtN)
# Delete "X" column
PresenceDatComb.df4 <- subset(PresenceDatComb.df4, select=-c(X))
#str(PresenceDatComb.df4)
## Group Descriptions into Habitat
PresenceDatComb.df4$Habitat <- PresenceDatComb.df4$Description
unique(PresenceDatComb.df4$Habitat)
PresenceDatComb.df4$Habitat[grepl("Rural ROW", PresenceDatComb.df4$Habitat)] <- "Rural Road ROW"
## Combine other categories
# Combine all Habitats with "Yard"
PresenceDatComb.df4$Habitat[grepl("Yard", PresenceDatComb.df4$Habitat)] <- "Yard"
# Combine all Habitats with "Trail ROW"
PresenceDatComb.df4$Habitat[grepl("Trail ROW", PresenceDatComb.df4$Habitat)] <- "Trail ROW"
# Combine all Habitats without "ROW" and "Yard"
PresenceDatComb.df4$Habitat[!grepl("ROW|Yard", PresenceDatComb.df4$Habitat)] <- "Field/Meadow"
# Combine Trail and Utility ROWs
PresenceDatComb.df4$Habitat[grepl("Trail|Utility", PresenceDatComb.df4$Habitat)] <- "Trail/Utility ROW"
# Combine Suburban and Urban ROWs
PresenceDatComb.df4$Habitat[grepl("Suburban|Urban", PresenceDatComb.df4$Habitat)] <- "Suburban/Urban Road ROW"
#
# Save data
write.table(PresenceDatComb.df4, sep = ",", col.names=NA, file=paste0(Species, "_Thinned", SpatFiltBuff, Units,"_presNLCD_", "1kPopDenIndexdat.csv"))
##################################################
## Aggregate by PopDenNHBRWtN based on Habitat and NewClassName
HabitatSummary.df <- aggregate(PopDenNHBRWtN ~ NewClassName + Habitat, data=PresenceDatComb.df4, sum)
sum(HabitatSummary.df$PopDenNHBRWtN)
########### Land Cover Class Between Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, Habitat ~ NewClassName, 
                   value.var = c("PopDenNHBRWtN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$Habitat
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassBetweenHabitatFrequencyTablePopDenNHBRWtN.csv"))
########### Land Cover Class Within Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, NewClassName ~ Habitat, 
                   value.var = c("PopDenNHBRWtN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$NewClassName
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassWithinHabitatFrequencyTablePopDenNHBRWtN.csv"))
##############################################################
## Aggregate by PopDenNHBRWtHN based on Habitat and NewClassName
HabitatSummary.df <- aggregate(PopDenNHBRWtHN ~ NewClassName + Habitat, data=PresenceDatComb.df4, sum)
sum(HabitatSummary.df$PopDenNHBRWtHN)
########### Land Cover Class Between Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, Habitat ~ NewClassName,
                   value.var = c("PopDenNHBRWtHN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$Habitat
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassBetweenHabitatFrequencyTablePopDenNHBRWtHN.csv"))
########### Land Cover Class Within Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, NewClassName ~ Habitat,
                   value.var = c("PopDenNHBRWtHN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$NewClassName
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassWithinHabitatFrequencyTablePopDenNHBRWtHN.csv"))
##############################################################
## Aggregate by RoadDistNHBRWtN based on Habitat and NewClassName
HabitatSummary.df <- aggregate(RoadDistNHBRWtN ~ NewClassName + Habitat, data=PresenceDatComb.df4, sum)
sum(HabitatSummary.df$RoadDistNHBRWtN)
########### Land Cover Class Between Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, Habitat ~ NewClassName,
                   value.var = c("RoadDistNHBRWtN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$Habitat
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassBetweenHabitatFrequencyTableRoadDistNHBRWtN.csv"))
########### Land Cover Class Within Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, NewClassName ~ Habitat,
                   value.var = c("RoadDistNHBRWtN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$NewClassName
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassWithinHabitatFrequencyTableRoadDistNHBRWtN.csv"))
##############################################################
## Aggregate by RoadDistNHBRWtHN based on Habitat and NewClassName
HabitatSummary.df <- aggregate(RoadDistNHBRWtHN ~ NewClassName + Habitat, data=PresenceDatComb.df4, sum)
sum(HabitatSummary.df$RoadDistNHBRWtHN)
########### Land Cover Class Between Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, Habitat ~ NewClassName,
                   value.var = c("RoadDistNHBRWtHN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$Habitat
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassBetweenHabitatFrequencyTableRoadDistNHBRWtHN.csv"))
########### Land Cover Class Within Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, NewClassName ~ Habitat,
                   value.var = c("RoadDistNHBRWtHN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$NewClassName
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassWithinHabitatFrequencyTableRoadDistNHBRWtHN.csv"))
##############################################################
## Aggregate by PopRoadHBRWtN based on Habitat and NewClassName
HabitatSummary.df <-aggregate(PopRoadHBRWtN ~ NewClassName + Habitat, data=PresenceDatComb.df4, sum)
sum(HabitatSummary.df$PopRoadHBRWtN)
########### Land Cover Class Between Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, Habitat ~ NewClassName,
                   value.var = c("PopRoadHBRWtN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$Habitat
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassBetweenHabitatFrequencyTablePopRoadHBRWtN.csv"))
########### Land Cover Class Within Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, NewClassName ~ Habitat,
                   value.var = c("PopRoadHBRWtN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$NewClassName
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassWithinHabitatFrequencyTablePopRoadHBRWtN.csv"))
##################
## Aggregate by PopRoadHBRWtN based on Habitat and NewClassName
HabitatSummary.df <-aggregate(PopRoadHBRWtHN ~ NewClassName + Habitat, data=PresenceDatComb.df4, sum)
sum(HabitatSummary.df$PopRoadHBRWtHN)
########### Land Cover Class Between Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, Habitat ~ NewClassName,
                   value.var = c("PopRoadHBRWtHN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$Habitat
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassBetweenHabitatFrequencyTablePopRoadHBRWtHN.csv"))
########### Land Cover Class Within Habitat Types
## Reshape from long to wide format
HabitatSummary.dfdcast <- dcast(HabitatSummary.df, NewClassName ~ Habitat,
                   value.var = c("PopRoadHBRWtHN"))
# Change NA values to zero
HabitatSummary.dfdcast[is.na(HabitatSummary.dfdcast)]=0
# Make first column as rownames
RowNames <- HabitatSummary.dfdcast$NewClassName
HabitatSummary.dfdcast <- HabitatSummary.dfdcast[,-1]
rownames(HabitatSummary.dfdcast) <- RowNames
## Make data frame of proportions
HabitatSummaryProp.df <- as.data.frame.matrix(prop.table(as.matrix(HabitatSummary.dfdcast), margin=2))
#
Sums.df1 <- data.frame(t(colSums(HabitatSummary.dfdcast)))
rowSums(Sums.df1)[1]
colnames(Sums.df1) <- colnames(HabitatSummaryProp.df)
rownames(Sums.df1)[1] <- "N"
HabitatSummaryProp.df2 <- rbind(HabitatSummaryProp.df, Sums.df1)
#str(HabitatSummaryProp.df2)
#sum(HabitatSummaryProp.df2[1:5,1])
# Save data with land cover classes
write.table(HabitatSummaryProp.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassWithinHabitatFrequencyTablePopRoadHBRWtHN.csv"))
##################
########### Land Cover Class Between Habitat Types
# Tabulate Raw Data counts on Habitat types based on land cover name
HabitatTable <- with(PresenceDatComb.df4, table(Habitat, NewClassName))
HabitatTable.df <- as.data.frame.matrix(addmargins(HabitatTable))
# Extract column sums for sample sizes
Sums.df <- data.frame(HabitatTable.df[(nrow(HabitatTable.df)),1:(ncol(HabitatTable.df)-1)])
rownames(Sums.df)[1] <- "N"
colnames(Sums.df) <- colnames(HabitatTable.df)[1:(ncol(HabitatTable.df)-1)]
## Make data frame of proportions
HabitatTableProp.df <- as.data.frame.matrix(prop.table(HabitatTable, margin=2))
# Add row sums for n
HabitatTableProb.df2 <- rbind(HabitatTableProp.df, Sums.df)
# Save data with land cover classes
write.table(HabitatTableProb.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassBetweenHabitatFrequencyTableRawData.csv"))
############
########### Land Cover Class Within Habitat Types
# Tabulate Raw Data counts on Habitat types based on land cover name
HabitatTable <- with(PresenceDatComb.df4, table(NewClassName, Habitat))
HabitatTable.df <- as.data.frame.matrix(addmargins(HabitatTable))
# Extract column sums for sample sizes
Sums.df <- data.frame(HabitatTable.df[(nrow(HabitatTable.df)),1:(ncol(HabitatTable.df)-1)])
rownames(Sums.df)[1] <- "N"
colnames(Sums.df) <- colnames(HabitatTable.df)[1:(ncol(HabitatTable.df)-1)]
## Make data frame of proportions
HabitatTableProp.df <- as.data.frame.matrix(prop.table(HabitatTable, margin=2))
# Add row sums for n
HabitatTableProb.df2 <- rbind(HabitatTableProp.df, Sums.df)
# Save data with land cover classes
write.table(HabitatTableProb.df2, sep = ",", col.names=NA, file=paste0(Species, "_LandCoverClassWithinHabitatFrequencyTableRawData.csv"))
####################################################################################
### Calculate accuracy assessment contingency table comparing Class Name with Revised Class Name
setwd(InDirect)
#
PresenceDatComb.df4 <- data.frame(read.csv(paste0(Species, "_Thinned", SpatFiltBuff, Units,"_presNLCD_", "1kPopDenIndexdat.csv"), stringsAsFactors=FALSE))
#
head(PresenceDatComb.df4)
#
RefPnts.mat <- as.matrix(PresenceDatComb.df4$RevClassName)
PredPnts.mat <- as.matrix(PresenceDatComb.df4$ClassName)
## Specify directory where function "ClassContingencyTable_Function.R" is kept
FunctDirectIn <- FunctDirect
Blank <- ""
RunType <- "Overall"
ModelNameOut <- ""
SubsetVariableNumber <- ""
VariableSubsets <- ""
output1 <- ""
Conditions2 <- ""
OutDirectIn3 <- InDirect
TestDataType <- ""
Run <- ""
#
ConfusionMat <- ClassContingencyTable2(RefPnts.mat, PredPnts.mat, Species, ModelNameOut, SubsetVariableNumber, VariableSubsets, output1,
Conditions2, OutDirectIn3, TestDataType, Run, RunType, FunctDirectIn, Blank, VarImp=FALSE, Output=TRUE)
#######################################################################################
# Calculate accuracy assessments by Habitat Type
HabitatList <-  unique(PresenceDatComb.df4$Habitat)
for(HabitatType in HabitatList) {
  #HabitatType <- HabitatList[5]
  PresenceDatComb.df4a <- PresenceDatComb.df4[which(PresenceDatComb.df4$Habitat==HabitatType),]
  head(PresenceDatComb.df4a)
  nrow(PresenceDatComb.df4a)
  RefPnts.mat <- as.matrix(PresenceDatComb.df4a$RevClassName)
  PredPnts.mat <- as.matrix(PresenceDatComb.df4a$ClassName)
  ## Specify directory where function "ClassContingencyTable_Function.R" is kept
  FunctDirectIn <- FunctDirect
  Blank <- ""
  RunType <- gsub("/", "", HabitatType)
  ModelNameOut <- ""
  SubsetVariableNumber <- ""
  VariableSubsets <- ""
  output1 <- ""
  Conditions2 <- ""
  OutDirectIn3 <- InDirect
  TestDataType <- ""
  Run <- ""
  #
  ConfusionMat <- ClassContingencyTable2(RefPnts.mat, PredPnts.mat, Species, ModelNameOut, SubsetVariableNumber, VariableSubsets, output1,
  Conditions2, OutDirectIn3, TestDataType, Run, RunType, FunctDirectIn, Blank, VarImp=FALSE, Output=TRUE)
}
