############################################### ############################################### ## GCWA Distance Analysis ## ## September 10, 2020 ## ############################################### ############################################### # This is an analysis of the 2018 point count # data on Golden-cheecked Warblers collected # at Balcones NWR using the gdistsamp function # in package 'unmarked' in program R. This analysis # uses layers with focal means at the 100 m scale. ####################################################### library(car) library(MuMIn) library(pls) library(Hmisc) library(ggcorrplot) library(unmarked) library(raster) ########################################################### ###### Read in and format data for gdistsamp ############## ########################################################### setwd("E:/GCWA") GCWA <- read.table("GCWA_RangewideDetection_ALL_20180614.csv", sep=",", header=TRUE) survey.info <- read.table("GCWA_RangewidePointCountSurveys_ALL_20190201.csv", sep=",", header=TRUE) survey.info$pt <- paste("PT.", survey.info$Point_number, sep="") wind <- tapply(survey.info$Wind_speed_code, list(survey.info$pt), mean) wind <- data.frame(wind=wind, pt=names(wind)) obs <- unique(survey.info[,c(5,35)]) obs2 <- unique(survey.info[,c(4,35)]) survey.data <- merge(wind, obs, by="pt") survey.data <- merge(survey.data, obs2, by="pt") # remove cleared points cleared <- c(23, 63, 99, 103, 783, 1233, 1256, 1507, 1671, 1792) GCWA <- GCWA[ !(GCWA$Point_number %in% cleared), ] cov <- read.table("cov_all_new_9_16_20.csv", sep=",", header=TRUE) cov.orig <- cov model_ht <- lm(cov$veg_ht_1000m_NAs_filled ~ cov$cc_3m_1000) cov$ht_diff <- resid(model_ht) cov.names <- names(cov) cov.names <- cov.names[!cov.names %in% c("Point_number", "point", "pt_ID")] for(i in 1:length(cov.names)){ cov[is.na(cov[,cov.names[i]]),cov.names[i]] <- median(cov[,cov.names[i]], na.rm = TRUE) } cov$Point_number <- factor(cov$Point_number) cov <- cov[ !(cov$Point_number %in% cleared), ] cov$jccper.100 <- cov$juniper_100m_pct.img/cov$treecov_100m_pct.img cov$jccper.1000 <- cov$juniper_1000m_pct.img/cov$treecov_1000m_pct.img cov$jccper.5000 <- cov$juniper_5000m_pct.img/cov$treecov_5000m_pct.img #### 2 m cov$jccper.100.2 <- (cov$jcc_2m_100)/cov$cc_2m_100 cov$jccper.1000.2 <- (cov$jcc_2m_1000)/cov$cc_2m_1000 cov$jccper.5000.2 <- (cov$jcc_2m_5000)/cov$cc_2m_5000 #### 3 m cov$jccper.100.3 <- (cov$jcc_3m_100)/cov$cc_3m_100 cov$jccper.1000.3 <- (cov$jcc_3m_1000)/cov$cc_3m_1000 cov$jccper.5000.3 <- (cov$jcc_3m_5000)/cov$cc_3m_5000 #### 4 m cov$jccper.100.4 <- (cov$jcc_4m_100)/cov$cc_4m_100 cov$jccper.1000.4 <- (cov$jcc_4m_1000)/cov$cc_4m_1000 cov$jccper.5000.4 <- (cov$jcc_4m_5000)/cov$cc_4m_5000 #### 5 m cov$jccper.100.5 <- (cov$jcc_5m_100)/cov$cc_5m_100 cov$jccper.1000.5 <- (cov$jcc_5m_1000)/cov$cc_5m_1000 cov$jccper.5000.5 <- (cov$jcc_5m_5000)/cov$cc_5m_5000 # 2m cov$cc_100_above2 <- cov$jcc_2m_100/cov$treecov_100m_pct.img cov$cc_1000_above2 <- cov$jcc_2m_1000/cov$treecov_1000m_pct.img cov$cc_5000_above2 <- cov$jcc_2m_5000/cov$treecov_5000m_pct.img # 3m cov$cc_100_above3 <- cov$jcc_3m_100/cov$treecov_100m_pct.img cov$cc_1000_above3 <- cov$jcc_3m_100/cov$treecov_1000m_pct.img cov$cc_5000_above3 <- cov$jcc_3m_100/cov$treecov_5000m_pct.img # 4m cov$cc_100_above4 <- cov$jcc_4m_100/cov$treecov_100m_pct.img cov$cc_1000_above4 <- cov$jcc_4m_100/cov$treecov_1000m_pct.img cov$cc_5000_above4 <- cov$jcc_4m_100/cov$treecov_5000m_pct.img # 5m cov$cc_100_above5 <- cov$jcc_5m_100/cov$treecov_100m_pct.img cov$cc_1000_above5 <- cov$jcc_5m_100/cov$treecov_1000m_pct.img cov$cc_5000_above5 <- cov$jcc_5m_100/cov$treecov_5000m_pct.img ################################################################### ### Standardize variables by substracting mean and dividing by sd # ################################################################### cov.m <- cov[ , -which(names(cov) %in% c("Point_number", "point", "pt_ID", "Group_Name"))] cov.m <- as.matrix(cov.m) cov.m.for.scale <- cov.m cov.m <- stdize(cov.m, center = TRUE, scale = TRUE) cov.m <-as.data.frame(cov.m) ################################################ # Examine correlations ######################### ################################################ res2 <- rcorr(as.matrix(cov.m)) res2 res2$P # ++++++++++++++++++++++++++++ # flattenCorrMatrix # ++++++++++++++++++++++++++++ # cormat : matrix of the correlation coefficients # pmat : matrix of the correlation p-values flattenCorrMatrix <- function(cormat, pmat) { ut <- upper.tri(cormat) data.frame( row = rownames(cormat)[row(cormat)[ut]], column = rownames(cormat)[col(cormat)[ut]], cor =(cormat)[ut], p = pmat[ut] ) } flat.cor<-flattenCorrMatrix(res2$r, res2$P) subset(flat.cor, cor > 0.7) index <- names(cov.orig) %in% c("treecov_100m_pct.img", "cc_3m_100", "cc_2m_100", "cc_4m_100", "cc_5m_100", "veg_ht100m_NAs_filled") mydata.cor <- cov.orig[,index] mydata.cor <- na.omit(mydata.cor) corr <- round(cor(mydata.cor), 1) ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE) ######################################## cov.m$point <- cov$point cov.m$point <- factor(cov.m$point) cov.m$Point_number <- cov$Point_number cov.m$region <- cov$Group_Name ################################################################### ##### Format point count data ##################################### ################################################################### #GCWA<-subset(GCWA, Detection_distance<225) GCWA <- subset(GCWA, Interval_0_3_min=="Y") # truncates data to first 3 min GCWA$Pt.id <- paste("PT.", GCWA$Point_number, sep="") bird.dist <- data.frame(point=GCWA$Pt.id, distance=GCWA$Detection_distance) bird.dist$distance <- as.numeric(as.character(bird.dist$distance)) # ensures distance is numeric bird.dist$point <- factor(bird.dist$point, levels=levels(as.factor(cov.m$point))) # adds in pts with no observations numbers as factors ################################################### ### Format data for unmarked ###################### ################################################### humps <- bird.dist[,2] %in% seq(0, 220, by=5) 151/length(humps) quantile(bird.dist[,2], 0.97) hist(bird.dist[,2], breaks=50) abline(v=25, col="red") abline(v=50, col="red") abline(v=75, col="red") cutpt <- c(0, 35, 66, 88, 110, 132, 154, 180)# Select cut points to minimize humping, max distance 189 m y2 <- formatDistData(bird.dist, distCol="distance", transectNameCol="point", dist.breaks=cutpt) y2 <- y2[order(rownames(y2)),] ########################################################### #### Create unmarked data frame ######################## ########################################################### cov.m <- cov.m[order(cov.m$point),] survey.data <- survey.data[order(survey.data$pt),] survey.data1 <- survey.data[!duplicated(survey.data$pt),] cov.m2 <- merge(cov.m, survey.data1, by.x="point", by.y="pt", all.x=TRUE, all.y=FALSE) cov.m2$wind[is.na(cov.m2$wind)] <- mean(cov.m2$wind, na.rm=TRUE) cov.m2$obs.simp2 <- factor(cov.m2$obs.simp2) cov.m2$obs.simp <- factor(cov.m2$obs.simp) umf <- unmarkedFrameGDS(y=y2, numPrimary=1, siteCovs=cov.m2, survey="point",dist.breaks=cutpt, unitsIn="m") summary(umf) #summary of compiled data hist(umf, xlab="distance (m)", main="", cex.lab=0.8, cex.axis=0.8) ############################################# ##Testing of different base-level functions # ############################################# basemodels <- list() basemodels$half.normal.nb <- gdistsamp(pformula=~obs.simp2+ruf100m_clp.img+treecov_100m_pct.img, phiformula=~1, lambdaformula= ~treecov_100m_pct.img+I(jccper.100^2)+I(srad_1000_clp.img^2), umf, output="density",rel.tol=0.00001, keyfun="halfnorm", mixture = "NB") basemodels$hazard.nb <- gdistsamp(pformula=~obs.simp2+ruf100m_clp.img+treecov_100m_pct.img, phiformula=~1, lambdaformula=~treecov_100m_pct.img+I(jccper.100^2)+I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.00001, keyfun="hazard", mixture = "NB") basemodels$half.normal.pois <- gdistsamp(pformula=~obs.simp2+ruf100m_clp.img+treecov_100m_pct.img, phiformula=~1, lambdaformula= ~treecov_100m_pct.img+I(jccper.100^2)+I(srad_1000_clp.img^2), umf, output="density",rel.tol=0.00001, keyfun="halfnorm") basemodels$hazard.pois <- gdistsamp(pformula=~obs.simp2+ruf100m_clp.img+treecov_100m_pct.img, phiformula=~1, lambdaformula=~treecov_100m_pct.img+I(jccper.100^2)+I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.00001, keyfun="hazard") ## Assemble the various model fits into a "fitList" and do model selection fits.basemodels <- fitList(fits=basemodels) (ms.basemodels <- modSel(fits.basemodels)) #nPars AIC delta AICwt cumltvWt #half.normal.pois 9 4852.25 0.00 0.7066 0.71 #half.normal.nb 10 4854.25 2.00 0.2599 0.97 #hazard.pois 10 4858.96 6.71 0.0246 0.99 #hazard.nb 11 4861.00 8.75 0.0089 1.00 ################################################### # No need for NB ################################## ################################################# # need add siteCovs=cov.m back in umf <- unmarkedFrameDS(y=y2, siteCovs=cov.m2, survey="point",dist.breaks=cutpt, unitsIn="m") summary(umf) #summary of compiled data hist(umf, xlab="distance (m)", main="", cex.lab=0.8, cex.axis=0.8) ############################################# ## Comparison of cc (orig, 2m, 3m, 4m, 5m) # ############################################# # Compare canopy cover, defined by height break, e.g. % cover of canopy at least 3 m cc.mods <- list() cc.mods$orig <- distsamp(~ obs.simp ~ treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods$cc2m <- distsamp(~obs.simp ~ cc_2m_100 + cc_2m_1000 + I(jccper.1000.2^2)+I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods$cc3m <- distsamp(~obs.simp ~ cc_3m_100 + cc_3m_1000 + I(jccper.1000.3^2)+I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods$cc4m <- distsamp(~obs.simp ~ cc_4m_100 + cc_4m_1000 + I(jccper.1000.4^2)+I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods$cc5m <- distsamp(~obs.simp ~ cc_5m_100 + cc_5m_1000 + I(jccper.1000.5^2)+I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods$null <- distsamp(~obs.simp ~ 1, data=umf, output="density", rel.tol=0.001) fits.cc.mods <- fitList(fits=cc.mods) (ms.cc.mods <- modSel(fits.cc.mods)) # Compare canopy cover, as % above given height break, e.g. of the total canopy cover 60% is above 3m cc.mods2 <- list() cc.mods2$orig <- distsamp(~obs.simp ~ treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods2$abovecc2m <- distsamp(~obs.simp ~ cc_100_above2 + treecov_100m_pct.img + cc_5000_above2 + treecov_1000m_pct.img + I(jccper.100.2) + I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods2$abovecc3m <- distsamp(~obs.simp ~ cc_100_above3 + treecov_100m_pct.img + cc_1000_above3 + treecov_1000m_pct.img + I(jccper.1000.3^2) + I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods2$abovecc4m <- distsamp(~obs.simp ~ cc_100_above4 + treecov_100m_pct.img + cc_1000_above4 + treecov_1000m_pct.img + I(jccper.1000.4^2) + I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) cc.mods2$abovecc5m <- distsamp(~obs.simp ~ cc_100_above5 + treecov_100m_pct.img + cc_1000_above5 + treecov_1000m_pct.img + I(jccper.1000.5^2) + I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) fits.cc.mods2 <- fitList(fits=cc.mods2) (ms.cc.mods2 <- modSel(fits.cc.mods2)) # % of canopy above 3m has most support. ############################################# ##Testing of different detection covariates # ############################################# detect.mods <- list() detect.mods$null <- distsamp(~1 ~treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), data=umf, output="density", rel.tol=0.001) detect.mods$high1km <- distsamp(~highways_1000_focal ~treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) detect.mods$allroads1km <- distsamp(~all_roads_1000_focal ~ treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) detect.mods$cc.100 <- distsamp(~treecov_100m_pct.img~ treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) detect.mods$cc.1000 <- distsamp(~treecov_1000m_pct.img ~ treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) detect.mods$ruf.100 <- distsamp(~ruf100m_clp.img ~ treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) detect.mods$ruf.1000 <- distsamp(~ruf1000m_clp.img ~ treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) detect.mods$wind <- distsamp(~wind ~treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) # detect.mods$obs.simp <- distsamp(~obs.simp ~treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) # detect.mods$ruf_all_roads1km <- distsamp(~ruf1000m_clp.img + all_roads_1000_focal ~treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) # detect.mods$ruf_allroads_obs.simp <- distsamp(~obs.simp + all_roads_1000_focal+ruf1000m_clp.img ~treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) # detect.mods$allroads_obs.simp <- distsamp(~obs.simp + all_roads_1000_focal ~treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) # detect.mods$ruf_obs.simp <-distsamp(~obs.simp + ruf1000m_clp.img ~treecov_100m_pct.img + treecov_1000m_pct.img + I(jccper.1000^2) + I(srad_1000_clp.img^2), umf, output="density", rel.tol=0.001) # fits.detect.mods <- fitList(fits=detect.mods) (ms.detect.mods <- modSel(fits.detect.mods)) ################################################## ##### Landscape level vegetation ########################### ################################################# veg.mods <- list() variab <- names(cov) ind1 <- grep("5km", variab) ind2 <- grep("jccper.5000", variab) ind3 <- grep("treecov_5000m_pct.img", variab) ind4 <- grep("_5000_", variab) index <- c(ind1, ind2, ind3, ind4) variab <- variab[index] variab <- variab[(c(1:9, 29, 39:50))] for(i in 1:length(variab)){ print(i) veg.mods[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ ", paste(variab[i]))), data=umf, output="density", rel.tol=0.001) # names(veg.mods)[i] <- paste(variab[i]) } for(i in (length(veg.mods)+1):(length(veg.mods)+length(variab))){ print(i) veg.mods[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ ", paste(variab[i-length(variab)], "+I(", variab[i-length(variab)], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(veg.mods)[i] <- paste(variab[i-length(variab)], "_", variab[i-length(variab)], "2", sep="") } for(i in (length(veg.mods)+1):(length(veg.mods)+length(variab))){ print(i) veg.mods[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ ", paste("treecov_5000m_pct.img", "*", variab[i-length(variab)*2], "+I(", variab[i-length(variab)*2], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(veg.mods)[i]<-paste("treecov_5000m_pct.img", "_x_", variab[i-length(variab)*2], "_", variab[i-length(variab)*2], "2", sep="") } for(i in (length(veg.mods)+1):(length(veg.mods)+length(variab))){ print(i) veg.mods[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ ", paste("treecov_5000m_pct.img", "*", variab[i-length(variab)*3]))), data=umf, output="density", rel.tol=0.001) # names(veg.mods)[i] <- paste("treecov_5000m_pct.img", "_x_", variab[i-length(variab)*3], sep="") } veg.mods$null <- distsamp(~obs.simp+all_roads_1000_focal+ruf1000m_clp.img ~1, data=umf, output="density", rel.tol=0.001) # fits.mod1 <- fitList(fits=veg.mods) (ms.mod1 <- modSel(fits.mod1)) veg.mods$treecov_5000m_pct.img_x_tin2017_5km.tif_tin2017_5km.tif2 <- NULL #Insig params treecov_5000m_pct.img_x_dur2017_5km.tif_dur2017_5km.tif2 <- NULL #insig params veg.mods$cc_5000_above3 <- distsamp(~obs.simp + ruf1000m_clp.img ~ treecov_5000m_pct.img + cc_5000_above3 + I(treecov_5000m_pct.img^2), data=umf, output="density", rel.tol=0.001) veg.mods$treecov_5000m_pct.img_x_eosn_ave_5km_cc_5000_above3_treecov_5000m_pct.img2_jccper.50002 <- distsamp(~obs.simp + ruf1000m_clp.img ~ treecov_5000m_pct.img*eosn_ave_5km + cc_5000_above3 + I(treecov_5000m_pct.img^2) , data=umf, output="density", rel.tol=0.001) fits.mod1 <- fitList(fits=veg.mods) (ms.mod1 <- modSel(fits.mod1)) summary(veg.mods$treecov_5000m_pct.img_x_eosn_ave_5km_cc_5000_above3_treecov_5000m_pct.img2_jccper.50002) ################################################## ##### Landscape level topo ########################### ################################################# topo.mods <- list() variab <- names(cov) ind1 <- grep("imi_100m_clp.img.1", variab) ind2 <- grep("5000_clp.img", variab) ind3 <- grep("5000m_clp.img", variab) index <- c(ind1, ind2, ind3) index <- unique(index) variab <- variab[index] variab <- variab[c(2:5, 7)] for(i in 1:length(variab)){ print(i) topo.mods[[i]] <- distsamp(formula=as.formula(paste("~obs.simp+all_roads_1000_focal+ruf1000m_clp.img ~ treecov_5000m_pct.img*eosn_ave_5km + cc_5000_above3 + I(treecov_5000m_pct.img^2) +", paste(variab[i]))), data=umf, output="density", rel.tol=0.001) # names(topo.mods)[i] <- paste("treecov_5000m_pct.img_x_cc_5000_above3", variab[i]) } for(i in 1:length(variab)){ print(i) topo.mods <- append(topo.mods, distsamp(formula=as.formula(paste("~obs.simp+all_roads_1000_focal+ruf1000m_clp.img ~ treecov_5000m_pct.img*eosn_ave_5km + cc_5000_above3 + I(treecov_5000m_pct.img^2) +", paste(variab[i], "+I(", variab[i], "^2)"))), data=umf, output="density", rel.tol=0.001)) # names(topo.mods)[length(topo.mods)] <- paste("treecov_5000m_pct.img_x_cc_5000_above3", variab[i], "_", variab[i], "2", sep="") } for(i in 1:length(variab)){ print(i) topo.mods <- append(topo.mods, distsamp(formula=as.formula(paste("~obs.simp+all_roads_1000_focal+ruf1000m_clp.img ~ treecov_5000m_pct.img*eosn_ave_5km + cc_5000_above3 + I(treecov_5000m_pct.img^2) +", paste("I(", variab[i], "^2)"))), data=umf, output="density", rel.tol=0.001)) # names(topo.mods)[length(topo.mods)] <- paste("treecov_5000m_pct.img_x_cc_5000_above3", variab[i], "2", sep="") } topo.mods$null <- distsamp(~ obs.simp + all_roads_1000_focal + ruf1000m_clp.img ~treecov_5000m_pct.img*cc_5000_above3, data=umf, output="density", rel.tol=0.001) # fits.mod1 <- fitList(fits=topo.mods) (ms.mod1 <- modSel(fits.mod1)) topo.mods$treecov_5000m_pct.img_x_cc_5000_above3el_5000m_clp.img_el_5000m_clp.img2 <- NULL #insig params topo.mods$treecov_5000m_pct.img_x_cc_5000_above3slope_deg_5000_clp.img_slope_deg_5000_clp.img2 <- NULL #insig params topo.mods$treecov_5000m_pct.img_x_cc_5000_above3cti_5000m_clp.img_cti_5000m_clp.img2 <- NULL #insig params fits.mod1 <- fitList(fits=topo.mods) (ms.mod1 <- modSel(fits.mod1)) ######################################################### ##### Landscape level climate ########################### ######################################################### climate.mods <- list() variab <- names(cov) ind1 <- grep("_5km", variab) index <- c(ind1) index <- unique(index) variab <- variab[index] variab <- variab[c(10:21, 25:26, 28)] for(i in 1:length(variab)){ print(i) climate.mods[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_5000m_pct.img*eosn_ave_5km + cc_5000_above3 + I(treecov_5000m_pct.img^2) + I(el_5000m_clp.img^2) + +", paste(variab[i]))), data=umf, output="density", rel.tol=0.001) # names(climate.mods)[i] <- paste(variab[i]) } for(i in (length(climate.mods)+1):(length(climate.mods)+length(variab))){ print(i) climate.mods[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_5000m_pct.img*eosn_ave_5km + cc_5000_above3 + I(treecov_5000m_pct.img^2) + I(el_5000m_clp.img^2) +", paste(variab[i-length(variab)], "+I(", variab[i-length(variab)], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(climate.mods)[i] <- paste(variab[i-length(variab)], "_", variab[i-length(variab)], "2", sep="") } for(i in (length(climate.mods)+1):(length(climate.mods)+length(variab))){ print(i) climate.mods[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_5000m_pct.img*eosn_ave_5km + cc_5000_above3 + I(treecov_5000m_pct.img^2) + I(el_5000m_clp.img^2) +", paste("I(", variab[i-length(variab)*2], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(climate.mods)[i] <- paste( variab[i-length(variab)*2], "2", sep="") } climate.mods$null <- distsamp(~obs.simp + ruf1000m_clp.img ~ treecov_5000m_pct.img*eosn_ave_5km + cc_5000_above3 + I(treecov_5000m_pct.img^2) + el_5000m_clp.img, data=umf, output="density", rel.tol=0.001) # fits.mod1 <- fitList(fits=climate.mods) (ms.mod1 <- modSel(fits.mod1)) climate.mods$temp_range_5km_temp_range_5km2 <- NULL climate.mods$temp_seasonality_5km_temp_seasonality_5km2 <- NULL fits.mod1 <- fitList(fits=climate.mods) (ms.mod1 <- modSel(fits.mod1)) ######################################### ### Site level veg ###################### ######################################### variab <- names(cov) ind1 <- grep("100m", variab) ind2 <- grep("1km", variab) ind3 <- grep("1000m", variab) ind4 <- grep("_1000_", variab) ind5 <- grep(".1000", variab) ind6 <- grep(".100", variab) index <- c(ind1, ind2, ind3, ind4, ind5, ind6) index <- unique(index) variab <- variab[index] variab <- variab[c(1, 7, 9, 17:24, 37, 40, 49:52, 53, 84:87)] mods.site.veg <- list() for(i in 1:length(variab)){ print(i) mods.site.veg[[i]] <- distsamp(formula=as.formula(paste("~obs.simp+ruf1000m_clp.img ~ ", paste(variab[i]))), data=umf, output="density", rel.tol=0.001) # names(mods.site.veg)[i] <- paste(variab[i]) } for(i in (length(mods.site.veg)+1):(length(mods.site.veg)+length(variab))){ print(i) mods.site.veg[[i]] <- distsamp(formula=as.formula(paste("~obs.simp+ruf1000m_clp.img ~ ", paste(variab[i-length(variab)], "+I(", variab[i-length(variab)], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(mods.site.veg)[i] <- paste(variab[i-length(variab)], "_", variab[i-length(variab)], "2", sep="") } for(i in (length(mods.site.veg)+1):(length(mods.site.veg)+length(variab))){ print(i) mods.site.veg[[i]] <- distsamp(formula=as.formula(paste("~obs.simp+ruf1000m_clp.img ~ ", paste("treecov_1000m_pct.img", "*", variab[i-length(variab)*2], "+I(", variab[i-length(variab)*2], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(mods.site.veg)[i]<-paste("treecov_1000m_pct.img", "_x_", variab[i-length(variab)*2], "_", variab[i-length(variab)*2], "2", sep="") } mods.site.veg$null <- distsamp(~obs.simp+ruf1000m_clp.img ~1, data=umf, output="density", rel.tol=0.001) # mods.site.veg$treecov_1000m_pct.img_x_eost_ave_1km_eost_ave_1km2 <- NULL #insig params fits.mod1 <- fitList(fits=mods.site.veg) (ms.mod1 <- modSel(fits.mod1)) ######################################### ### Site level topo ########################## ######################################### variab <- names(cov) ind1 <- grep("100m", variab) ind2 <- grep("1km", variab) ind3 <- grep("1000m", variab) ind4 <- grep("_1000_", variab) ind5 <- grep(".1000", variab) ind6 <- grep(".100", variab) index <- c(ind1, ind2, ind3, ind4, ind5, ind6) index <- unique(index) variab <- variab[index] variab <- variab[c(2, 4, 6, 8, 13, 14, 33:34, 36, 43, 44, 45, 46, 68)] mods.site.topo <- list() for(i in 1:length(variab)){ print(i) mods.site.topo[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img + treecov_100m_pct.img+cc_100_above3 + I(jccper.1000^2)+", paste(variab[i]))), data=umf, output="density", rel.tol=0.001) # names(mods.site.topo)[i] <- paste(variab[i]) } for(i in (length(mods.site.topo)+1):(length(mods.site.topo)+length(variab))){ print(i) mods.site.topo[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img + treecov_100m_pct.img+cc_100_above3 + I(jccper.1000^2)+ ", paste(variab[i-length(variab)], "+I(", variab[i-length(variab)], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(mods.site.topo)[i] <- paste(variab[i-length(variab)], "_", variab[i-length(variab)], "2", sep="") } for(i in (length(mods.site.topo)+1):(length(mods.site.topo)+length(variab))){ print(i) mods.site.topo[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img + treecov_100m_pct.img+cc_100_above3 + I(jccper.1000^2)+ ", paste("I(", variab[i-length(variab)*2], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(mods.site.topo)[i] <- paste(variab[i-length(variab)*2], "2", sep="") } for(i in (length(mods.site.topo)+1):(length(mods.site.topo)+length(variab))){ print(i) mods.site.topo[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img + treecov_100m_pct.img+cc_100_above3 + I(jccper.1000^2)+ treecov_1000m_pct.img*", variab[i-length(variab)*3])), data=umf, output="density", rel.tol=0.001) # names(mods.site.topo)[i] <- paste("treecov_1000m_pct.img_x_",variab[i-length(variab)*3], sep="") } mods.site.topo$null <- distsamp(~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img + treecov_100m_pct.img+cc_100_above3 + I(jccper.1000^2), data=umf, output="density", rel.tol=0.001) # mods.site.topo$el_1000m_clp.img_el_1000m_clp.img2 <- NULL #insig params mods.site.topo$intel_100m_clp.img <-NULL #insig params mods.site.topo$cc_1000_above3_cc_1000_above32 <-NULL #insig params mods.site.topo$treecov_1000m_pct.img_x_el_100m_clp.img <-NULL #insig params mods.site.topo$treecov_1000m_pct.img_x_el_1000m_clp.img <-NULL #insig params mods.site.topo$treecov_1000m_pct.img_x_slope_deg_1000_clp.img <-NULL #insig params mods.site.topo$cc_1000_above2 <-NULL #insig params mods.site.topo$el_100m_clp.img_el_100m_clp.img2 <- NULL #insig params mods.site.topo$cc_1000_above2_cc_1000_above22 <- NULL #insig params mods.site.topo$cc_1000_above3 <- NULL #insig params mods.site.topo$cti_1000m_clp.img_cti_1000m_clp.img2 <- NULL #insig params fits.mod1 <- fitList(fits=mods.site.topo) (ms.mod1 <- modSel(fits.mod1)) ######################################### ### Site/territory level anthro ######### ######################################### variab <- names(cov) ind1 <- grep("100m", variab) ind2 <- grep("1km", variab) ind3 <- grep("1000m", variab) ind4 <- grep("_1000_", variab) ind5 <- grep(".1000", variab) ind6 <- grep(".100", variab) index <- c(ind1, ind2, ind3, ind4, ind5, ind6) index <- unique(index) variab <- variab[index] variab <- variab[c(5, 10:12, 35, 41:42, 46:47, 69:70)] mods.site.anthro <- list() for(i in 1:length(variab)){ print(i) mods.site.anthro[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img *cti_1000m_clp.img + treecov_100m_pct.img+cc_100_above3 + I(jccper.1000^2)+", paste(variab[i]))), data=umf, output="density", rel.tol=0.001) # names(mods.site.anthro)[i] <- paste(variab[i]) } for(i in (length(mods.site.anthro)+1):(length(mods.site.anthro)+length(variab))){ print(i) mods.site.anthro[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img *cti_1000m_clp.img + treecov_100m_pct.img+cc_100_above3 + I(jccper.1000^2)+", paste(variab[i-length(variab)], "+I(", variab[i-length(variab)], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(mods.site.anthro)[i] <- paste(variab[i-length(variab)], "_", variab[i-length(variab)], "2", sep="") } for(i in (length(mods.site.anthro)+1):(length(mods.site.anthro)+length(variab))){ print(i) mods.site.anthro[[i]] <- distsamp(formula=as.formula(paste("~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img *cti_1000m_clp.img + treecov_100m_pct.img+cc_100_above3 + I(jccper.1000^2) +", paste("I(", variab[i-length(variab)*2], "^2)"))), data=umf, output="density", rel.tol=0.001) # names(mods.site.anthro)[i] <- paste(variab[i-length(variab)*2], "2", sep="") } mods.site.anthro$null <- distsamp(~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img * cti_1000m_clp.img + treecov_100m_pct.img + cc_100_above3 + I(jccper.1000^2) , data=umf, output="density", rel.tol=0.001) # mods.site.anthro$all_roads_1000_focal_all_roads_1000_focal2 <- NULL mods.site.anthro$forloss_ld_100m_clp.img_forloss_ld_100m_clp.img2 <- NULL mods.site.anthro$forloss_dens_100mcir_clp.tif_forloss_dens_100mcir_clp.tif2 <- NULL mods.site.anthro$srad_1000_clp.img_srad_1000_clp.img2 <- NULL mods.site.anthro$build_dens_100m_clp.tif2 <- NULL mods.site.anthro$forloss_ld_1000m_clp.img2 <- NULL mods.site.anthro$build_dens_1000m_clp.tif <- NULL mods.site.anthro$build_dens_100m_clp.tif_build_dens_100m_clp.tif2 <- NULL mods.site.anthro$forloss_dens_100mcir_clp.tif <- NULL mods.site.anthro$srad_1000_clp.img <- NULL mods.site.anthro$build_dis_100m_clp.tif <- NULL mods.site.anthro$forloss_ld_1000m_clp.img <- NULL mods.site.anthro$all_roads_1000_focal <- NULL mods.site.anthro$all_roads_1000_focal2 <- NULL mods.site.anthro$build_dens_1000m_clp.tif2 <- NULL mods.site.anthro$highways_100_focal <- NULL mods.site.anthro$build_dis_1000m_clp.tif <- NULL mods.site.anthro$all_roads2_100_focal <- NULL mods.site.anthro$build_dens_1000m_clp.tif_build_dens_1000m_clp.tif2 <- NULL mods.site.anthro$forloss_ld_1000m_clp.img_forloss_ld_1000m_clp.img2 <- NULL mods.site.anthro$highways_100_focal_highways_100_focal2 <- NULL mods.site.anthro$all_roads2_100_focal_all_roads2_100_focal2 <- NULL mods.site.anthro$highways_100_focal2 <- NULL mods.site.anthro$all_roads2_100_focal2 <- NULL mods.site.anthro$srad_1000_clp.img2 <- NULL mods.site.anthro$forloss_ld_100m_clp.img <- NULL # duplicate of forloss mods.site.anthro$forloss_dens_100mcir_clp.tif2 <- NULL # duplicate of forloss mods.site.anthro$build_dis_1000m_clp.tif_build_dis_1000m_clp.tif2 <- NULL # Insig params mods.site.anthro$build_dis_100m_clp.tif_build_dis_100m_clp.tif2 <- NULL # Insig params mods.site.anthro$build_dis_1000m_clp.tif_build_dis_1000m_clp.tif2 <- NULL mods.site.anthro$build_dis_100m_clp.tif_build_dis_100m_clp.tif2 <- NULL fits.mod1 <- fitList(fits=mods.site.anthro) (ms.mod1 <- modSel(fits.mod1)) ############################################################ # Variables from earlier steps that had significant support # ############################################################ variables <- data.frame(cov$treecov_1000m_pct.img, cov$temp_range_5km, cov$cti_1000m_clp.img, cov$el_100m_clp.img, cov$treecov_100m_pct.img, cov$cc_100_above3, cov$jccper.1000.3^2, cov$build_dis_1000m_clp.tif^2, cov$forloss_ld_100m_clp.img^2, cov$srad_1000_clp.img^2, cov$treecov_5000m_pct.img, cov$eosn_ave_5km, cov$el_5000m_clp.img, cov$el_5000m_clp.img^2, cov$el_1000m_clp.img^2) ######################################################### ## Examine correlations among variables ################# ######################################################### # Computing correlation matrix correlation_matrix <- round(cor(variables),2) # Plotting correlation matrix ggcorrplot(correlation_matrix, hc.order = TRUE, type = "lower", lab = TRUE) global <- distsamp(formula = ~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img + I(temp_range_5km^2) + cti_1000m_clp.img * treecov_1000m_pct.img + treecov_100m_pct.img + cc_100_above3 + I(jccper.1000.3^2) + I(build_dis_1000m_clp.tif^2) + I(srad_1000_clp.img^2) + treecov_5000m_pct.img*eosn_ave_5km + I(el_5000m_clp.img^2) + I(el_1000m_clp.img^2), data = umf, output = "density", rel.tol = 0.001) ### All possible subsets ### Warning - the below is very slow, to see final model skip to final model below dd <- dredge(global) top.models <- get.models(dd, cumsum(weight) <= 0.99) # inspect models, remove those with non-sig terms (>0.15) and/or those with correlated (>0.7) terms index <- 1 for(i in 1:length(top.models)){ print(i) if (all(summary(top.models[[i]])$coefficients[,4][-1] <0.15)){ index <- c(index, i) } else { print(i) } } index <- index[-1] length(index) final_model <- distsamp(formula = ~obs.simp + ruf1000m_clp.img ~ treecov_1000m_pct.img + I(temp_range_5km^2) + cti_1000m_clp.img*treecov_1000m_pct.img + treecov_100m_pct.img + cc_100_above3 + I(jccper.1000.3^2) , data = umf, output = "density", rel.tol = 0.001) summary(final_model) ############################################ ## Create spatial predictions ############## ############################################ cov.m.for.scale <- as.data.frame(cov.m.for.scale) setwd("G:/Layers_used_for_prediction") jcc <- raster("jccper_1000_scaled_re.tif") build <- raster("buildden_1000_scaled_re_c_aoi_final.tif") temp <- raster("temp_scaled_c_aoi_final.tif") cti <- raster("cti_scaled_c_aoi_final.tif") cc_100 <- raster("cc_100_scaled_re_c_aoi_final3.tif") cc_1000 <- raster("cc_1000_scaled_re_c_aoi_final3.tif") sost <- raster("sost_scaled_c_aoi_final.tif") cc_per <- raster("cc_100_above3_scaled_c_aoi_final2.tif") srad <- raster("srad_1000_scaled_re_c_aoi_final2.tif") p_stack <- stack(cti, temp, cc_1000, cc_100, jcc, cc_per) names(p_stack) <- c("cti_1000m_clp.img", "temp_range_5km", "treecov_1000m_pct.img", "treecov_100m_pct.img", "jccper.1000.3", "cc_100_above3") model.predict <- function(x){ predict(final_model, type="state", newdata=x, na.action = na.pass) } xseq <- seq(334908.2, 682328.2, length.out = 100) #Use foreach loop and %dopar% command for(i in 1:99){ print(i) alli <- crop(p_stack, extent(xseq[i], xseq[i+1], 3234670, 3662050)) print("crop done") p <- data.frame(rasterToPoints(alli)) print("raster to points done") names(p) summary(p) if(dim(p)[1]==0)next chunk1 <- seq(1,length(p[,1]), by=1000) startTime=date() result.1 <- foreach(a=chunk1)%dopar% { cat('\r',a/dim(p)[1]*100); flush.console() samp.data <- p[a:(a+999),] samp.data <- as.data.frame(samp.data) samp.data <- samp.data[!is.na(samp.data[,1]),] samp.data2 <- samp.data samp.data2[is.na(samp.data2)] <- 0 data_out <- model.predict(samp.data2) samp.data$den <- data_out[,1] samp.data$den[rowSums(is.na(samp.data[,3:9]))>0] <- NA samp.data$lower <- data_out[,3] samp.data$lower[rowSums(is.na(samp.data[,3:9]))>0] <- NA samp.data$upper <- data_out[,4] samp.data$upper[rowSums(is.na(samp.data[,3:9]))>0] <- NA samp.data3 <- samp.data[,c("x", "y", "den", "lower", "upper")] return(samp.data3) } gcw.mod <- do.call(rbind, result.1) # Make list into data frame colnames(gcw.mod) <- c("x", "y", "den", "lower", "upper") gcw.mod_den <- rasterFromXYZ(gcw.mod[,c("x", "y", "den")]) # Make raster of predicted density gcw.mod_low <- rasterFromXYZ(gcw.mod[,c("x", "y", "lower")]) # Make raster of predicted density gcw.mod_high <- rasterFromXYZ(gcw.mod[,c("x", "y", "upper")]) # Make raster of predicted density setwd("G:/new_output") writeRaster(gcw.mod_den, filename=paste("gcwa_prediction_no_sost_den_3min", i, ".tif", sep="_"),format="GTiff",datatype="FLT4S", overwrite=TRUE) writeRaster(gcw.mod_low, filename=paste("gcwa_prediction_no_sost_lower_3min", i, ".tif", sep="_"),format="GTiff",datatype="FLT4S", overwrite=TRUE) writeRaster(gcw.mod_high, filename=paste("gcwa_prediction_no_sost_upperz_3min", i, ".tif", sep="_"),format="GTiff",datatype="FLT4S", overwrite=TRUE) }