From 6a1378eb46f279d21c14d55802df7a15d0a61e0d Mon Sep 17 00:00:00 2001 From: Guilherme Gall Date: Sat, 1 Apr 2017 22:19:51 -0300 Subject: [PATCH] Reformata os arquivos em ENM/fct/ (issue #1) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Usei o formatR para colocar os códigos mais próximos do estilo recomendado pelo Wickham. A função usada foi tidy_dir('ENM/fct/', width.cutoff=80, indent=2) Comentários inline e definição de funções foram alteradas manualmente. --- ENM/fct/createBuffer.R | 80 +- ENM/fct/cropModel.R | 10 +- ENM/fct/eigenvariables.fct.R | 57 +- ENM/fct/ensemble.R | 139 +-- ENM/fct/final_model.R | 245 +++--- ENM/fct/modelos.R | 1572 +++++++++++++++++----------------- ENM/fct/read.eval1.R | 351 ++++---- ENM/fct/write_Raster.R | 13 +- 8 files changed, 1221 insertions(+), 1246 deletions(-) diff --git a/ENM/fct/createBuffer.R b/ENM/fct/createBuffer.R index b8c993c..43c2d4b 100644 --- a/ENM/fct/createBuffer.R +++ b/ENM/fct/createBuffer.R @@ -1,43 +1,39 @@ -createBuffer <- function(coord_ = coord, - sp_ = sp, - occs_ = occs, - seed_ = seed, - n.back_ = n.back, - buffer.type_ = "mean" - ){ -#Transformando em spatial points -coordinates(coord_) <- ~lon+lat - -if (buffer.type_ == "mean") dist.buf <- mean(spDists(x = coord_, longlat = FALSE, segments = TRUE)) -if (buffer.type_ == "max") dist.buf <- max(spDists(x = coord_, longlat = FALSE, segments = TRUE)) - -buffer <- raster::buffer(coord_, width = dist.buf, dissolve = TRUE) - -#Transformando coords de novo em matriz para rodar resto script -coord_ <- occs_[occs_$sp == sp_,c('lon','lat')] - -#Transformando em spatial polygon data frame -buffer <- SpatialPolygonsDataFrame(buffer,data=as.data.frame(buffer@plotOrder), match.ID = FALSE) -crs(buffer) <- crs(predictors) - -#########TENHO CERTEZA DE QUE ISTO PODE FICAR MENOS PESADO -#Reference raster com mesmo extent e resolution que predictors -r_buffer <- crop(predictors, buffer) -r_buffer <- mask(r_buffer, buffer) - -#r_buffer <- raster(ext=extent(predictors), resolution=res(predictors)) - -#Rasterizando o buffer p/ geração dos ptos aleatorios -#r_buffer <- rasterize(buffer, r_buffer, field=buffer@plotOrder) -#Limitando a mascara ambiental -#r_buffer <- r_buffer*(predictors[[1]]!=0) -################ - - -#Gerando pontos aleatorios no buffer -set.seed(seed_+2) -backgr <- randomPoints(r_buffer, n.back_) -rm(buffer) -gc() -return(backgr) +createBuffer <- function(coord_ = coord, sp_ = sp, occs_ = occs, seed_ = seed, n.back_ = n.back, + buffer.type_ = "mean") { + # Transformando em spatial points + coordinates(coord_) <- ~lon + lat + + if (buffer.type_ == "mean") + dist.buf <- mean(spDists(x = coord_, longlat = FALSE, segments = TRUE)) + if (buffer.type_ == "max") + dist.buf <- max(spDists(x = coord_, longlat = FALSE, segments = TRUE)) + + buffer <- raster::buffer(coord_, width = dist.buf, dissolve = TRUE) + + # Transformando coords de novo em matriz para rodar resto script + coord_ <- occs_[occs_$sp == sp_, c("lon", "lat")] + + # Transformando em spatial polygon data frame + buffer <- SpatialPolygonsDataFrame(buffer, data = as.data.frame(buffer@plotOrder), + match.ID = FALSE) + crs(buffer) <- crs(predictors) + + ######### TENHO CERTEZA DE QUE ISTO PODE FICAR MENOS PESADO Reference raster com mesmo + ######### extent e resolution que predictors + r_buffer <- crop(predictors, buffer) + r_buffer <- mask(r_buffer, buffer) + + # r_buffer <- raster(ext=extent(predictors), resolution=res(predictors)) + + # Rasterizando o buffer p/ geração dos ptos aleatorios r_buffer <- + # rasterize(buffer, r_buffer, field=buffer@plotOrder) Limitando a mascara + # ambiental r_buffer <- r_buffer*(predictors[[1]]!=0) + + + # Gerando pontos aleatorios no buffer + set.seed(seed_ + 2) + backgr <- randomPoints(r_buffer, n.back_) + rm(buffer) + gc() + return(backgr) } diff --git a/ENM/fct/cropModel.R b/ENM/fct/cropModel.R index f3ba192..edeef63 100644 --- a/ENM/fct/cropModel.R +++ b/ENM/fct/cropModel.R @@ -1,5 +1,5 @@ -cropModel <- function(modelo, mascara){ - modelo <- crop(modelo, mascara) - modelo <- mask(modelo, mascara) -return(modelo) - } \ No newline at end of file +cropModel <- function(modelo, mascara) { + modelo <- crop(modelo, mascara) + modelo <- mask(modelo, mascara) + return(modelo) +} diff --git a/ENM/fct/eigenvariables.fct.R b/ENM/fct/eigenvariables.fct.R index c23029f..27bf3d3 100644 --- a/ENM/fct/eigenvariables.fct.R +++ b/ENM/fct/eigenvariables.fct.R @@ -1,42 +1,33 @@ -eigenvariables.fct <- function(vars, name, proportion = 0.95){ +eigenvariables.fct <- function(vars, name, proportion = 0.95) { library("raster") - if (file.exists("./pca") == FALSE) dir.create("./pca") - #Running PCA: - # Counts not NA cels + if (file.exists("./pca") == FALSE) + dir.create("./pca") + # Running PCA: Counts not NA cels non.na <- sum(!is.na(values(vars[[1]]))) # Sample the study area with n-non.na and creates an environmental table sr <- sampleRandom(vars, non.na) # faz o PCA dessa tabela padronizada pca <- prcomp(scale(sr)) summary.pca <- summary(pca) - #Saving results: - capture.output(pca, file = sprintf('./pca/%s.pca.txt', name)) - - #saving summary - capture.output(summary.pca, file = sprintf('./pca/%s.summary.pca.txt', name)) - - #Plotting results - #GGPLOT - ##### - #library(ggplot2) - # create data frame with scores - #scores <- as.data.frame(pca$x) - # plot of observations - #ggplot(data = scores, aes(x = PC1, y = PC2, label = rownames(scores))) + - # geom_hline(yintercept = 0, colour = "gray65") + - # geom_vline(xintercept = 0, colour = "gray65") + - # geom_text(colour = "tomato", alpha = 0.8, size = 4) + - # ggtitle("PCA plot of USA States - Crime Rates") - ##### - # png(filename = sprintf('./pca/%s.pca.biplot.png',name), - # bg = "white") - # biplot(pca) - # dev.off() - + # Saving results: + capture.output(pca, file = sprintf("./pca/%s.pca.txt", name)) + + # saving summary + capture.output(summary.pca, file = sprintf("./pca/%s.summary.pca.txt", name)) + + # Plotting results GGPLOT library(ggplot2) create data frame with scores scores + # <- as.data.frame(pca$x) plot of observations ggplot(data = scores, aes(x = PC1, + # y = PC2, label = rownames(scores))) + geom_hline(yintercept = 0, colour = + # 'gray65') + geom_vline(xintercept = 0, colour = 'gray65') + geom_text(colour = + # 'tomato', alpha = 0.8, size = 4) + ggtitle('PCA plot of USA States - Crime + # Rates') png(filename = sprintf('./pca/%s.pca.biplot.png',name), bg = 'white') + # biplot(pca) dev.off() + # Creating eigenvariable in space - axis.nb <- which(summary.pca$importance["Cumulative Proportion",] >= proportion)[1] + axis.nb <- which(summary.pca$importance["Cumulative Proportion", ] >= proportion)[1] eigenvariables <- predict(vars, pca, index = 1:axis.nb) - - if (file.exists("./env") == FALSE) dir.create("./env") - writeRaster(eigenvariables,sprintf('./env/%s.eigenvariables.tif',name),overwrite=T) -} \ No newline at end of file + + if (file.exists("./env") == FALSE) + dir.create("./env") + writeRaster(eigenvariables, sprintf("./env/%s.eigenvariables.tif", name), overwrite = T) +} diff --git a/ENM/fct/ensemble.R b/ENM/fct/ensemble.R index f914042..c7bba47 100644 --- a/ENM/fct/ensemble.R +++ b/ENM/fct/ensemble.R @@ -1,71 +1,78 @@ -ensemble <- function(sp, - models.dir = "./models",##onde estao os modelos em geral ("models") - final.dir = "presfinal", ##onde estao os modelos finais ("presfinal") - ensemble.dir = "ensemble", - occs = spp.filt, - which.models = c("Final.bin.mean3","Final.mean.bin7"), - consensus = F, - consensus.level = 0.5#cuanto de los modelos sea retenido: 0.5majority - ){ +## onde estao os modelos em geral ('models') onde estao os modelos finais +## ('presfinal') cuanto de los modelos sea retenido: 0.5majority - ##pasta de output - if (file.exists(paste0(models.dir,"/",sp,"/",ensemble.dir,"/")) == FALSE){ - dir.create(paste0(models.dir,"/",sp,"/",ensemble.dir,"/")) - }# - library(raster) - library(scales) - library(maps) - - ## para cada tipo de modelo - for (whi in which.models){ - - cat(paste(whi,"-",sp,"\n")) #lê os arquivos - tif.files <- list.files(paste0(models.dir,"/",sp,"/",final.dir),full.names=T,pattern=paste0(whi,'.*tif$')) - - if(length(tif.files)==0) { - cat(paste("No models to ensemble from for",sp,"\n")) - } else { - cat(paste(length(tif.files),"models to ensemble from for",sp,"\n")) - mod2 <- raster::stack(tif.files) - if(length(tif.files)==1){ - ensemble.m <- mod2 - } else { - #plot(mod2) +ensemble <- function(sp, models.dir = "./models", final.dir = "presfinal", ensemble.dir = "ensemble", + occs = spp.filt, which.models = c("Final.bin.mean3", "Final.mean.bin7"), consensus = F, + consensus.level = 0.5) { + + ## pasta de output + if (file.exists(paste0(models.dir, "/", sp, "/", ensemble.dir, "/")) == FALSE) + { + dir.create(paste0(models.dir, "/", sp, "/", ensemble.dir, "/")) + } # + library(raster) + library(scales) + library(maps) + + ## para cada tipo de modelo + for (whi in which.models) { + + cat(paste(whi, "-", sp, "\n")) #lê os arquivos + tif.files <- list.files(paste0(models.dir, "/", sp, "/", final.dir), full.names = T, + pattern = paste0(whi, ".*tif$")) + + if (length(tif.files) == 0) { + cat(paste("No models to ensemble from for", sp, "\n")) + } else { + cat(paste(length(tif.files), "models to ensemble from for", sp, "\n")) + mod2 <- raster::stack(tif.files) + if (length(tif.files) == 1) { + ensemble.m <- mod2 + } else { + # plot(mod2) ensemble.m <- mean(mod2) - ensemble.sd <- raster::overlay(mod2,fun=function(x){return(sd(x,na.rm=T))}) - } - coord <- occs[occs$sp==sp,c('lon','lat')] - - #par(mar=c(5,4,1,1)) - #plot(ensemble.m,main=paste(sp,whi,"ensemble"),font.main=3) - #points(coord,pch=19,col=alpha("grey60",0.6)) - #map('world',c('',"South America"),add=T,col="grey") - - png(filename=paste0(models.dir,"/",sp,"/",ensemble.dir,"/",sp,"_",whi,"_ensemble.png"),res=300,width=410*300/72,height=480*300/72) - par(mfrow=c(1,1),mar=c(3,4,4,0)) - plot(ensemble.m,main=paste(sp,whi),legend=F, - cex.main=1,font.main=3) - map('world',c('',"South America"),add=T,col="grey") - points(coord,pch=21,cex=0.6,bg=scales::alpha("cyan",0.6)) + ensemble.sd <- raster::overlay(mod2, fun = function(x) { + return(sd(x, na.rm = T)) + }) + } + coord <- occs[occs$sp == sp, c("lon", "lat")] + + # par(mar=c(5,4,1,1)) plot(ensemble.m,main=paste(sp,whi,'ensemble'),font.main=3) + # points(coord,pch=19,col=alpha('grey60',0.6)) map('world',c('','South + # America'),add=T,col='grey') + + png(filename = paste0(models.dir, "/", sp, "/", ensemble.dir, "/", sp, + "_", whi, "_ensemble.png"), res = 300, width = 410 * 300/72, height = 480 * + 300/72) + par(mfrow = c(1, 1), mar = c(3, 4, 4, 0)) + plot(ensemble.m, main = paste(sp, whi), legend = F, cex.main = 1, font.main = 3) + map("world", c("", "South America"), add = T, col = "grey") + points(coord, pch = 21, cex = 0.6, bg = scales::alpha("cyan", 0.6)) + dev.off() + + # o ensemble cru + writeRaster(ensemble.m, filename = paste0(models.dir, "/", sp, "/", ensemble.dir, + "/", sp, "_", whi, "_ensemble.tif"), overwrite = T) + + #### Consensus models + if (consensus == TRUE) { + ensemble.consensus <- ensemble.m >= consensus.level + writeRaster(ensemble.consensus, filename = paste0(models.dir, "/", + sp, "/", ensemble.dir, "/", sp, "_", whi, "_ensemble", consensus.level * + 100, ".tif"), overwrite = T) + + + png(filename = paste0(models.dir, "/", sp, "/", ensemble.dir, "/", + sp, "_", whi, "_ensemble", consensus.level * 100, ".png"), res = 300, + width = 410 * 300/72, height = 480 * 300/72) + par(mfrow = c(1, 1), mar = c(3, 4, 4, 0)) + plot(ensemble.consensus, main = paste(whi, consensus.level * 100), + legend = F, cex.main = 1, font.main = 3) + map("world", c("", "South America"), add = T, col = "grey") + points(coord, pch = 19, cex = 0.3, col = scales::alpha("cyan", 0.6)) dev.off() - - # o ensemble cru - writeRaster(ensemble.m,filename=paste0(models.dir,"/",sp,"/",ensemble.dir,"/",sp,"_",whi,"_ensemble.tif"),overwrite=T) - - ####Consensus models - if (consensus == TRUE){ - ensemble.consensus <- ensemble.m >= consensus.level - writeRaster(ensemble.consensus,filename=paste0(models.dir,"/",sp,"/",ensemble.dir,"/",sp,"_",whi,"_ensemble",consensus.level*100,".tif"),overwrite=T) - - - png(filename=paste0(models.dir,"/",sp,"/",ensemble.dir,"/",sp,"_",whi,"_ensemble",consensus.level*100,".png"),res=300,width=410*300/72,height=480*300/72) - par(mfrow=c(1,1),mar=c(3,4,4,0)) - plot(ensemble.consensus,main=paste(whi,consensus.level*100),legend=F,cex.main=1,font.main=3) - map('world',c('',"South America"),add=T,col="grey") - points(coord,pch=19,cex=0.3,col=scales::alpha("cyan",0.6)) - dev.off() - } + } } -} - #return(ensemble.m) + } + # return(ensemble.m) } diff --git a/ENM/fct/final_model.R b/ENM/fct/final_model.R index aa68d8c..52d4803 100644 --- a/ENM/fct/final_model.R +++ b/ENM/fct/final_model.R @@ -1,140 +1,113 @@ ######### Final modeling: one model per algorithm ---- -finalModel <- function(sp, - select.partitions = T, - algoritmos = c("maxent", "rf", "svm"),#NULL, - threshold = c("spec_sens"), - TSS.value = 0.7, - models.dir = "./models", - final.dir = "presfinal") { - if (file.exists(paste0(models.dir,"/",sp,"/",final.dir)) == FALSE) - dir.create(paste0(models.dir,"/",sp,"/",final.dir), recursive = TRUE) - print(date()) - - cat(paste(sp,"\n")) - library("raster") - library("data.table") - cat(paste("Reading the evaluation files","\n")) - evall <- - list.files(path = paste0(models.dir,"/",sp), pattern = "evaluate",full.names = T) - lista <- list() - for (i in 1:length(evall)) { - lista[[i]] <- read.table(file = evall[i],header = T,row.names = 1) - } - stats <- rbindlist(lista) - stats <- as.data.frame(stats) - - #Extracts only for the selected algorithm - if (exists("algoritmos")==F) algoritmos <- unique(stats$algoritmo) - algoritmos <- as.factor(algoritmos) - #algo <- algoritmos[1] - todo <- stack() - for (algo in algoritmos) { - cat(paste("Extracting data for",algo,"\n")) - stats2 <- stats[stats$algoritmo == algo,] - - - part <- nrow(stats2)#How many partitions were there - - cat(paste("Reading models from .tif files","\n")) - modelos.cont <- - list.files( - path = paste0(models.dir,"/",sp),full.names = T,pattern = paste0(algo,"_cont_") - ) - - modelos.bin <- - list.files( - path = paste0(models.dir,"/",sp),full.names = T,pattern = paste0(algo,"_bin_") - ) - mod.cont <- stack(modelos.cont)#(0) - - mod.bin <- stack(modelos.bin)#(0) - - names(mod.cont) <- paste0(sp,algo,"Partition",1:part) - - names(mod.bin) <- names(mod.cont) - - if (select.partitions == T) { - cat("selecting partitions for", sp, algo,"\n") - sel.index <- which(stats2[,"TSS"] >= TSS.value) - cont.sel <- mod.cont[[sel.index]]#(1) - bin.sel <- mod.bin[[sel.index]]#(5) - - - th.mean <- mean(stats2[,names(stats2) == threshold][sel.index]) - - if (length(sel.index) == 0) - cat(paste("No partition was selected for",sp,algo,"\n")) - - #en caso de que sea solo uno varios modelos son el mismo - if (length(sel.index) == 1) { - cat(paste( - length(sel.index), "partition was selected for",sp,algo,"\n" - )) - - final <- stack(bin.sel,#[3], - bin.sel#[7] - ) - names(final) <- - c( - "Final.bin.mean3", - "Final.mean.bin7" - ) - } - - #en caso de que sean más aplica el mapa - - if (length(sel.index) > 1) { - cat(paste( - length(sel.index), "partitions were selected for",sp,"\n" - )) - - final.cont.mean <- mean(cont.sel)#(2) - final.bin.mean <- (final.cont.mean > th.mean)#(3) - final.sel.bin <- mean(bin.sel)#(7) - - - final <- - stack(final.bin.mean, - final.sel.bin - ) - names(final) <- - c( - "Final.bin.mean3", - "Final.mean.bin7" - ) - } - - - if (exists("final")) { - #plot(final) - #Escribe final - rasterOptions(setfileext = F) - writeRaster( - x = final,filename = paste0( - models.dir,"/",sp,"/",final.dir,"/",names(final),sp,algo,".tif" - ),bylayer = T,overwrite = T,format = "GTiff" - ) - - for (i in 1:dim(final)[[3]]) { - png( - filename = paste0( - models.dir,"/",sp,"/",final.dir,"/",names(final)[i],sp,algo,".png" - ) - ) - plot(final[[i]],main = names(final)[i]) - dev.off() - } - todo <- addLayer(todo,final) - } - cat("select",sp,algo,"DONE","\n") - - } +finalModel <- function(sp, select.partitions = T, algoritmos = c("maxent", "rf", + "svm"), threshold = c("spec_sens"), TSS.value = 0.7, models.dir = "./models", + final.dir = "presfinal") { + if (file.exists(paste0(models.dir, "/", sp, "/", final.dir)) == FALSE) + dir.create(paste0(models.dir, "/", sp, "/", final.dir), recursive = TRUE) + print(date()) + + cat(paste(sp, "\n")) + library("raster") + library("data.table") + cat(paste("Reading the evaluation files", "\n")) + evall <- list.files(path = paste0(models.dir, "/", sp), pattern = "evaluate", + full.names = T) + lista <- list() + for (i in 1:length(evall)) { + lista[[i]] <- read.table(file = evall[i], header = T, row.names = 1) + } + stats <- rbindlist(lista) + stats <- as.data.frame(stats) + + # Extracts only for the selected algorithm + if (exists("algoritmos") == F) + algoritmos <- unique(stats$algoritmo) + algoritmos <- as.factor(algoritmos) + # algo <- algoritmos[1] + todo <- stack() + for (algo in algoritmos) { + cat(paste("Extracting data for", algo, "\n")) + stats2 <- stats[stats$algoritmo == algo, ] + + + part <- nrow(stats2) #How many partitions were there + + cat(paste("Reading models from .tif files", "\n")) + modelos.cont <- list.files(path = paste0(models.dir, "/", sp), full.names = T, + pattern = paste0(algo, "_cont_")) + + modelos.bin <- list.files(path = paste0(models.dir, "/", sp), full.names = T, + pattern = paste0(algo, "_bin_")) + mod.cont <- stack(modelos.cont) #(0) + + mod.bin <- stack(modelos.bin) #(0) + + names(mod.cont) <- paste0(sp, algo, "Partition", 1:part) + + names(mod.bin) <- names(mod.cont) + + if (select.partitions == T) { + cat("selecting partitions for", sp, algo, "\n") + sel.index <- which(stats2[, "TSS"] >= TSS.value) + cont.sel <- mod.cont[[sel.index]] #(1) + bin.sel <- mod.bin[[sel.index]] #(5) + + + th.mean <- mean(stats2[, names(stats2) == threshold][sel.index]) + + if (length(sel.index) == 0) + cat(paste("No partition was selected for", sp, algo, "\n")) + + # en caso de que sea solo uno varios modelos son el mismo + if (length(sel.index) == 1) { + cat(paste(length(sel.index), "partition was selected for", sp, algo, + "\n")) + + # bin.sel #[3] bin.sel #[7] + final <- stack(bin.sel, bin.sel) + names(final) <- c("Final.bin.mean3", "Final.mean.bin7") + } + + # en caso de que sean más aplica el mapa + + if (length(sel.index) > 1) { + cat(paste(length(sel.index), "partitions were selected for", sp, + "\n")) + + final.cont.mean <- mean(cont.sel) #(2) + final.bin.mean <- (final.cont.mean > th.mean) #(3) + final.sel.bin <- mean(bin.sel) #(7) + + + final <- stack(final.bin.mean, final.sel.bin) + names(final) <- c("Final.bin.mean3", "Final.mean.bin7") + } + + + if (exists("final")) { + # plot(final) Escribe final + rasterOptions(setfileext = F) + writeRaster(x = final, filename = paste0(models.dir, "/", sp, "/", + final.dir, "/", names(final), sp, algo, ".tif"), bylayer = T, overwrite = T, + format = "GTiff") + + for (i in 1:dim(final)[[3]]) { + png(filename = paste0(models.dir, "/", sp, "/", final.dir, "/", + names(final)[i], sp, algo, ".png")) + plot(final[[i]], main = names(final)[i]) + dev.off() } - - print(paste("DONE",algo,"\n")) - print(date()) - - - print(date()) + todo <- addLayer(todo, final) + } + cat("select", sp, algo, "DONE", "\n") + + } + } + + print(paste("DONE", algo, "\n")) + print(date()) + + + print(date()) } diff --git a/ENM/fct/modelos.R b/ENM/fct/modelos.R index 2911505..0340318 100644 --- a/ENM/fct/modelos.R +++ b/ENM/fct/modelos.R @@ -1,846 +1,860 @@ # Cada função de modelagem do_ALGORITMO retorna o threshold e escreve os arquivos # de saída no diretório informado pelo parâmetro models.dir -do_bioclim <- function(sp, - predictors = predictors, - pres_train = pres_train, #NEW - pres_test = pres_test, #NEW - backg_test = backg_test, #NEW - i = i, #NEW - models.dir = model.dir, - project.model = project.model, - projections = projections, - projdata = projdata,#um vector con nombres - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ){ - cat(paste("Bioclim",'\n')) - bc <- bioclim (predictors, pres_train) - ebc <- evaluate(pres_test,backg_test,bc,predictors) - thresholdbc <- ebc@t[which.max(ebc@TPR + ebc@TNR)] - thbc <- threshold(ebc) - bc_TSS <- max(ebc@TPR + ebc@TNR)-1 - - bc_cont <- predict(predictors, bc, progress='text') - bc_bin <- bc_cont > thresholdbc - bc_cut <- bc_cont * bc_bin - thbc$AUC <- ebc@auc - thbc$TSS <- bc_TSS - thbc$algoritmo <- "BioClim" - thbc$partition <- i - row.names(thbc) <- paste(sp,i,"BioClim") +# mask → a SpatialPolygonsDataFrame layer to mask and crop the predicted model - if (class(mask) == "SpatialPolygonsDataFrame"){ +do_bioclim <- function(sp, + predictors = predictors, + pres_train = pres_train, + pres_test = pres_test, + backg_test = backg_test, + i = i, + models.dir = model.dir, + project.model = project.model, + projections = projections, + projdata = projdata, + mask = mask) { + cat(paste("Bioclim", "\n")) + bc <- bioclim(predictors, pres_train) + ebc <- evaluate(pres_test, backg_test, bc, predictors) + thresholdbc <- ebc@t[which.max(ebc@TPR + ebc@TNR)] + thbc <- threshold(ebc) + bc_TSS <- max(ebc@TPR + ebc@TNR) - 1 + + bc_cont <- predict(predictors, bc, progress = "text") + bc_bin <- bc_cont > thresholdbc + bc_cut <- bc_cont * bc_bin + thbc$AUC <- ebc@auc + thbc$TSS <- bc_TSS + thbc$algoritmo <- "BioClim" + thbc$partition <- i + row.names(thbc) <- paste(sp, i, "BioClim") + + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + bc_cont <- cropModel(bc_cont, mask) + bc_bin <- cropModel(bc_bin, mask) + bc_cut <- cropModel(bc_cut, mask) + + } + writeRaster(x = bc_cont, filename = paste0(models.dir, "/", sp, "/BioClim_cont_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = bc_bin, filename = paste0(models.dir, "/", sp, "/BioClim_bin_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = bc_cut, filename = paste0(models.dir, "/", sp, "/BioClim_cut_", + sp, "_", i, ".tif"), overwrite = T) + + png(filename = paste0(models.dir, "/", sp, "/Bioclim", sp, "_", i, "%03d.png")) + plot(bc_cont, main = paste("Bioclim raw", "\n", "AUC =", round(ebc@auc, 2), "-", + "TSS =", round(bc_TSS, 2))) + plot(bc_bin, main = paste("Bioclim P/A", "\n", "AUC =", round(ebc@auc, 2), "-", + "TSS =", round(bc_TSS, 2))) + plot(bc_cut, main = paste("Bioclim cut", "\n", "AUC =", round(ebc@auc, 2), "-", + "TSS =", round(bc_TSS, 2))) + dev.off() + + if (project.model == T) { + for (proj in projections) { + # data <- list.files(paste0('./env/',proj),pattern=proj) data2 <- stack(data) + stopifnot(names(projdata) == names(predictors)) + bc_proj <- predict(projdata, bc, progress = "text") + bc_proj_bin <- bc_proj > thresholdbc + bc_proj_cut <- bc_proj_bin * bc_proj + # Normaliza o modelo cut bc_proj_cut <- bc_proj_cut/maxValue(bc_proj_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { source("./fct/cropModel.R") - bc_cont <- cropModel(bc_cont , mask) - bc_bin <- cropModel(bc_bin , mask) - bc_cut <- cropModel(bc_cut , mask) - + bc_proj <- cropModel(bc_proj, mask) + bc_proj_bin <- cropModel(bc_proj_bin, mask) + bc_proj_cut <- cropModel(bc_proj_cut, mask) + } + writeRaster(x = bc_proj, filename = paste0(models.dir, "/", sp, "/", + proj, "/BioClim_cont_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = bc_proj_bin, filename = paste0(models.dir, "/", sp, "/", + proj, "/BioClim_bin_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = bc_proj_cut, filename = paste0(models.dir, "/", sp, "/", + proj, "/BioClim_cut_", sp, "_", i, ".tif"), overwrite = T) + rm(data2) + png(filename = paste0(models.dir, "/", sp, "/", proj, "/Bioclim", sp, + "_", i, "%03d.png")) + plot(bc_proj, main = paste("Bioclim raw", "\n", "AUC =", round(ebc@auc, + 2), "-", "TSS =", round(bc_TSS, 2))) + plot(bc_proj_bin, main = paste("Bioclim P/A", "\n", "AUC =", round(ebc@auc, + 2), "-", "TSS =", round(bc_TSS, 2))) + plot(bc_proj_cut, main = paste("Bioclim cut", "\n", "AUC =", round(ebc@auc, + 2), "-", "TSS =", round(bc_TSS, 2))) + dev.off() } - writeRaster(x=bc_cont,filename=paste0(models.dir,"/",sp,"/BioClim_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=bc_bin,filename=paste0(models.dir,"/",sp,"/BioClim_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=bc_cut,filename=paste0(models.dir,"/",sp,"/BioClim_cut_",sp,"_",i,".tif"),overwrite=T) - - png(filename=paste0(models.dir,"/",sp,"/Bioclim",sp,"_",i,"%03d.png")) - plot(bc_cont,main=paste("Bioclim raw","\n","AUC =", round(ebc@auc,2),'-',"TSS =",round(bc_TSS,2))) - plot(bc_bin,main=paste("Bioclim P/A","\n","AUC =", round(ebc@auc,2),'-',"TSS =",round(bc_TSS,2))) - plot(bc_cut,main=paste("Bioclim cut","\n","AUC =", round(ebc@auc,2),'-',"TSS =",round(bc_TSS,2))) - dev.off() - - if (project.model == T){ - for (proj in projections){ - #data <- list.files(paste0("./env/",proj),pattern=proj) - #data2 <- stack(data) - stopifnot(names(projdata) == names(predictors)) - bc_proj <- predict(projdata,bc,progress='text') - bc_proj_bin <- bc_proj > thresholdbc - bc_proj_cut <- bc_proj_bin * bc_proj - # Normaliza o modelo cut - #bc_proj_cut <- bc_proj_cut/maxValue(bc_proj_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - bc_proj <- cropModel(bc_proj , mask) - bc_proj_bin <- cropModel(bc_proj_bin , mask) - bc_proj_cut <- cropModel(bc_proj_cut , mask) - } - writeRaster(x=bc_proj,filename=paste0(models.dir,"/",sp,"/",proj,"/BioClim_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=bc_proj_bin,filename=paste0(models.dir,"/",sp,"/",proj,"/BioClim_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=bc_proj_cut,filename=paste0(models.dir,"/",sp,"/",proj,"/BioClim_cut_",sp,"_",i,".tif"),overwrite=T) - rm(data2) - png(filename=paste0(models.dir,"/",sp,"/",proj,"/Bioclim",sp,"_",i,"%03d.png")) - plot(bc_proj,main=paste("Bioclim raw","\n","AUC =", round(ebc@auc,2),'-',"TSS =",round(bc_TSS,2))) - plot(bc_proj_bin,main=paste("Bioclim P/A","\n","AUC =", round(ebc@auc,2),'-',"TSS =",round(bc_TSS,2))) - plot(bc_proj_cut,main=paste("Bioclim cut","\n","AUC =", round(ebc@auc,2),'-',"TSS =",round(bc_TSS,2))) - dev.off() - } - } - return(thbc) + } + return(thbc) } do_randomForest <- function(sp, - predictors = predictors, - sdmdata_train = sdmdata_train, - envtest_pre = envtest_pre, - envtest_back = envtest_back, - i = i, - envtrain = envtrain, - models.dir = model.dir, - project.model = project.model, - projections = projections, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ){ - library(randomForest) - cat(paste("Random Forests",'\n')) - #rf1 <- tuneRF(x=envtrain,y=sdmdata_train$pa,stepFactor = 0.5) - rf <- randomForest (sdmdata_train$pa~.,data=envtrain) - #rf <- randomForest (x =envtrain ,y=factor(sdmdata_train$pa),xtest=envtest,ytest = factor(sdmdata_teste$pa))#fazendo teste interno a funcao evaluate nao serve :( - - erf <- evaluate(envtest_pre,envtest_back,rf) - rf_TSS <- max(erf@TPR + erf@TNR)-1 - - thresholdrf <- erf@t[which.max(erf@TPR + erf@TNR)] - thrf <- threshold(erf) - thrf$AUC <- erf@auc - thrf$TSS <- rf_TSS#raro - thrf$algoritmo <- "rf" - thrf$partition <- i - row.names(thrf) <- paste(sp,i,"rf") - - rf_cont <- predict(predictors,rf,progress='text',type="response") - rf_bin <- rf_cont>thresholdrf - rf_cut <- rf_bin * rf_cont - #rf1_cut <- rf1_cut/maxValue(rf1_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ + predictors = predictors, + sdmdata_train = sdmdata_train, + envtest_pre = envtest_pre, + envtest_back = envtest_back, + i = i, + envtrain = envtrain, + models.dir = model.dir, + project.model = project.model, + projections = projections, + mask = mask) { + library(randomForest) + cat(paste("Random Forests", "\n")) + # rf1 <- tuneRF(x=envtrain,y=sdmdata_train$pa,stepFactor = 0.5) + rf <- randomForest(sdmdata_train$pa ~ ., data = envtrain) + # rf <- randomForest (x =envtrain ,y=factor(sdmdata_train$pa),xtest=envtest,ytest + # = factor(sdmdata_teste$pa))#fazendo teste interno a funcao evaluate nao serve + # :( + + erf <- evaluate(envtest_pre, envtest_back, rf) + rf_TSS <- max(erf@TPR + erf@TNR) - 1 + + thresholdrf <- erf@t[which.max(erf@TPR + erf@TNR)] + thrf <- threshold(erf) + thrf$AUC <- erf@auc + thrf$TSS <- rf_TSS #raro + thrf$algoritmo <- "rf" + thrf$partition <- i + row.names(thrf) <- paste(sp, i, "rf") + + rf_cont <- predict(predictors, rf, progress = "text", type = "response") + rf_bin <- rf_cont > thresholdrf + rf_cut <- rf_bin * rf_cont + # rf1_cut <- rf1_cut/maxValue(rf1_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + rf_cont <- cropModel(rf_cont, mask) + rf_bin <- cropModel(rf_bin, mask) + rf_cut <- cropModel(rf_cut, mask) + } + writeRaster(x = rf_cont, filename = paste0(models.dir, "/", sp, "/rf_cont_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = rf_bin, filename = paste0(models.dir, "/", sp, "/rf_bin_", sp, + "_", i, ".tif"), overwrite = T) + writeRaster(x = rf_cut, filename = paste0(models.dir, "/", sp, "/rf_cut_", sp, + "_", i, ".tif"), overwrite = T) + + png(filename = paste0(models.dir, "/", sp, "/rf", sp, "_", i, "%03d.png")) + plot(rf_cont, main = paste("RF raw", "\n", "AUC =", round(erf@auc, 2), "-", "TSS =", + round(rf_TSS, 2))) + plot(rf_bin, main = paste("RF P/A", "\n", "AUC =", round(erf@auc, 2), "-", "TSS =", + round(rf_TSS, 2))) + plot(rf_cut, main = paste("RF cut", "\n", "AUC =", round(erf@auc, 2), "-", "TSS =", + round(rf_TSS, 2))) + dev.off() + + + if (project.model == T) { + for (proj in projections) { + data <- list.files(paste0("./env/", proj), pattern = proj) + data2 <- stack(data) + rf_proj <- predict(data2, rf, progress = "text") + rf_proj_bin <- rf_proj > thresholdrf + rf_proj_cut <- rf_proj_bin * rf_proj + # Normaliza o modelo cut rf_proj_cut <- rf_proj_cut/maxValue(rf_proj_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { source("./fct/cropModel.R") - rf_cont <- cropModel(rf_cont , mask) - rf_bin <- cropModel(rf_bin , mask) - rf_cut <- cropModel(rf_cut , mask) - } - writeRaster(x=rf_cont,filename=paste0(models.dir,"/",sp,"/rf_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=rf_bin,filename=paste0(models.dir,"/",sp,"/rf_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=rf_cut,filename=paste0(models.dir,"/",sp,"/rf_cut_",sp,"_",i,".tif"),overwrite=T) - - png(filename=paste0(models.dir,"/",sp,"/rf",sp,"_",i,"%03d.png")) - plot(rf_cont,main=paste("RF raw","\n","AUC =", round(erf@auc,2),'-',"TSS =",round(rf_TSS,2))) - plot(rf_bin,main=paste("RF P/A","\n","AUC =", round(erf@auc,2),'-',"TSS =",round(rf_TSS,2))) - plot(rf_cut,main=paste("RF cut","\n","AUC =", round(erf@auc,2),'-',"TSS =",round(rf_TSS,2))) - dev.off() - - - if (project.model == T){ - for (proj in projections){ - data <- list.files(paste0("./env/",proj),pattern=proj) - data2 <- stack(data) - rf_proj <- predict(data2,rf,progress='text') - rf_proj_bin <- rf_proj > thresholdrf - rf_proj_cut <- rf_proj_bin * rf_proj - # Normaliza o modelo cut - #rf_proj_cut <- rf_proj_cut/maxValue(rf_proj_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - rf_proj <- cropModel(rf_proj , mask) - rf_proj_bin <- cropModel(rf_proj_bin , mask) - rf_proj_cut <- cropModel(rf_proj_cut , mask) - } - writeRaster(x=rf_proj,filename=paste0(models.dir,"/",sp,"/",proj,"/rf_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=rf_proj_bin,filename=paste0(models.dir,"/",sp,"/",proj,"/rf_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=rf_proj_cut,filename=paste0(models.dir,"/",sp,"/",proj,"/rf_cut_",sp,"_",i,".tif"),overwrite=T) - rm(data2) - } + rf_proj <- cropModel(rf_proj, mask) + rf_proj_bin <- cropModel(rf_proj_bin, mask) + rf_proj_cut <- cropModel(rf_proj_cut, mask) + } + writeRaster(x = rf_proj, filename = paste0(models.dir, "/", sp, "/", + proj, "/rf_cont_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = rf_proj_bin, filename = paste0(models.dir, "/", sp, "/", + proj, "/rf_bin_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = rf_proj_cut, filename = paste0(models.dir, "/", sp, "/", + proj, "/rf_cut_", sp, "_", i, ".tif"), overwrite = T) + rm(data2) } - return(thrf) + } + return(thrf) } do_SVM <- function(sp, - predictors = predictors, - sdmdata_train = sdmdata_train, - envtest_pre = envtest_pre, - envtest_back = envtest_back, - i = i, - envtrain = envtrain, - part = part, - models.dir = model.dir, - project.model = project.model, - projections = projections, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ){ - cat(paste("SVM",'\n')) - library(kernlab) - svm <- ksvm(sdmdata_train$pa~.,data=envtrain,cross=part)##svm deve ser com a variável resposta binária ou contínua, eu acho que binária - esvm <- evaluate(envtest_pre,envtest_back,svm) - #esvm <- evaluate(pres_test,backg_test,model = svm,x = predictors) - svm_TSS <- max(esvm@TPR + esvm@TNR)-1 - thresholdsvm <- esvm@t[which.max(esvm@TPR + esvm@TNR)] - thsvm <- threshold (esvm) - thsvm$AUC <- esvm@auc - thsvm$TSS <- svm_TSS - thsvm$algoritmo <- "svm" - thsvm$partition <- i - row.names(thsvm) <- paste(sp,i,"svm") - svm_cont <- predict(predictors,svm,progress='text') - svm_bin <- svm_cont>thresholdsvm - svm_cut <- svm_bin * svm_cont - - #TRANSFORMA 0 A 1 - svm_cont <- svm_cont/maxValue(svm_cont) - svm_cut <- svm_cut/maxValue(svm_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ + predictors = predictors, + sdmdata_train = sdmdata_train, + envtest_pre = envtest_pre, + envtest_back = envtest_back, + i = i, + envtrain = envtrain, + part = part, + models.dir = model.dir, + project.model = project.model, + projections = projections, + mask = mask) { + cat(paste("SVM", "\n")) + library(kernlab) + svm <- ksvm(sdmdata_train$pa ~ ., data = envtrain, cross = part) ##svm deve ser com a variável resposta binária ou contínua, eu acho que binária + esvm <- evaluate(envtest_pre, envtest_back, svm) + # esvm <- evaluate(pres_test,backg_test,model = svm,x = predictors) + svm_TSS <- max(esvm@TPR + esvm@TNR) - 1 + thresholdsvm <- esvm@t[which.max(esvm@TPR + esvm@TNR)] + thsvm <- threshold(esvm) + thsvm$AUC <- esvm@auc + thsvm$TSS <- svm_TSS + thsvm$algoritmo <- "svm" + thsvm$partition <- i + row.names(thsvm) <- paste(sp, i, "svm") + svm_cont <- predict(predictors, svm, progress = "text") + svm_bin <- svm_cont > thresholdsvm + svm_cut <- svm_bin * svm_cont + + # TRANSFORMA 0 A 1 + svm_cont <- svm_cont/maxValue(svm_cont) + svm_cut <- svm_cut/maxValue(svm_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + svm_cont <- cropModel(svm_cont, mask) + svm_bin <- cropModel(svm_bin, mask) + svm_cut <- cropModel(svm_cut, mask) + } + + writeRaster(x = svm_cont, filename = paste0(models.dir, "/", sp, "/svm_cont_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = svm_bin, filename = paste0(models.dir, "/", sp, "/svm_bin_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = svm_cut, filename = paste0(models.dir, "/", sp, "/svm_cut_", + sp, "_", i, ".tif"), overwrite = T) + + png(filename = paste0(models.dir, "/", sp, "/svm", sp, "_", i, "%03d.png")) + plot(svm_cont, main = paste("SVM raw", "\n", "AUC =", round(esvm@auc, 2), "-", + "TSS =", round(svm_TSS, 2))) + plot(svm_bin, main = paste("SVM P/A", "\n", "AUC =", round(esvm@auc, 2), "-", + "TSS =", round(svm_TSS, 2))) + plot(svm_cut, main = paste("SVM cut", "\n", "AUC =", round(esvm@auc, 2), "-", + "TSS =", round(svm_TSS, 2))) + dev.off() + + if (project.model == T) { + for (proj in projections) { + data <- list.files(paste0("./env/", proj), pattern = proj) + data2 <- stack(data) + svm_proj <- predict(data2, svm, progress = "text") + svm_proj_bin <- svm_proj > thresholdsvm + svm_proj_cut <- svm_proj_bin * svm_proj + + # Normaliza o modelo cut svm_proj_cut <- svm_proj_cut/maxValue(svm_proj_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { source("./fct/cropModel.R") - svm_cont <- cropModel(svm_cont , mask) - svm_bin <- cropModel(svm_bin , mask) - svm_cut <- cropModel(svm_cut , mask) + svm_proj <- cropModel(svm_proj, mask) + svm_proj_bin <- cropModel(svm_proj_bin, mask) + svm_proj_cut <- cropModel(svm_proj_cut, mask) + } + writeRaster(x = svm_proj, filename = paste0(models.dir, "/", sp, "/", + proj, "/svm_cont_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = svm_proj_bin, filename = paste0(models.dir, "/", sp, + "/", proj, "/svm_bin_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = svm_proj_cut, filename = paste0(models.dir, "/", sp, + "/", proj, "/svm_cut_", sp, "_", i, ".tif"), overwrite = T) + rm(data2) } - - writeRaster(x=svm_cont,filename=paste0(models.dir,"/",sp,"/svm_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=svm_bin,filename=paste0(models.dir,"/",sp,"/svm_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=svm_cut,filename=paste0(models.dir,"/",sp,"/svm_cut_",sp,"_",i,".tif"),overwrite=T) - - png(filename=paste0(models.dir,"/",sp,"/svm",sp,"_",i,"%03d.png")) - plot(svm_cont,main=paste("SVM raw","\n","AUC =", round(esvm@auc,2),'-',"TSS =",round(svm_TSS,2))) - plot(svm_bin,main=paste("SVM P/A","\n","AUC =", round(esvm@auc,2),'-',"TSS =",round(svm_TSS,2))) - plot(svm_cut,main=paste("SVM cut","\n","AUC =", round(esvm@auc,2),'-',"TSS =",round(svm_TSS,2))) - dev.off() - - if (project.model == T){ - for (proj in projections){ - data <- list.files(paste0("./env/",proj),pattern=proj) - data2 <- stack(data) - svm_proj <- predict(data2,svm,progress='text') - svm_proj_bin <- svm_proj > thresholdsvm - svm_proj_cut <- svm_proj_bin * svm_proj - - # Normaliza o modelo cut - #svm_proj_cut <- svm_proj_cut/maxValue(svm_proj_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - svm_proj <- cropModel(svm_proj , mask) - svm_proj_bin <- cropModel(svm_proj_bin , mask) - svm_proj_cut <- cropModel(svm_proj_cut , mask) - } - writeRaster(x=svm_proj,filename=paste0(models.dir,"/",sp,"/",proj,"/svm_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=svm_proj_bin,filename=paste0(models.dir,"/",sp,"/",proj,"/svm_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=svm_proj_cut,filename=paste0(models.dir,"/",sp,"/",proj,"/svm_cut_",sp,"_",i,".tif"),overwrite=T) - rm(data2) - } - } - return(thsvm) + } + return(thsvm) } do_maxent <- function(sp, - predictors = predictors, - pres_train = pres_train, #NEW - pres_test = pres_test, #NEW - backg_test = backg_test, #NEW - i = i, #NEW - models.dir = model.dir, - project.model = project.model, - projections = projections, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ){ - cat(paste("maxent",'\n')) - Sys.setenv(NOAWT=TRUE)#descomentei para ver - library(rJava) - mx <- maxent (predictors, pres_train) - png(filename = paste0(models.dir,"/",sp,"/maxent_variable_contribution_",sp,"_",i,".png")) - plot(mx) - dev.off() - - png(filename = paste0(models.dir,"/",sp,"/maxent_variable_response_",sp,"_",i,".png")) - response(mx) - dev.off() - - emx <- evaluate(pres_test,backg_test,mx,predictors) - thresholdmx <- emx@t[which.max(emx@TPR + emx@TNR)] - thmx <- threshold(emx) - mx_TSS <- max(emx@TPR + emx@TNR)-1 - mx_cont <- predict(mx, predictors,progress='text') - mx_bin <- mx_cont > thresholdmx - mx_cut <- mx_cont * mx_bin - thmx$AUC <- emx@auc - thmx$TSS <- mx_TSS - thmx$algoritmo <- "maxent" - thmx$partition <- i - row.names(thmx) <- paste(sp,i,"maxent") - if (class(mask) == "SpatialPolygonsDataFrame"){ + predictors = predictors, + pres_train = pres_train, + pres_test = pres_test, + backg_test = backg_test, + i = i, + models.dir = model.dir, + project.model = project.model, + projections = projections, + mask = mask) { + cat(paste("maxent", "\n")) + Sys.setenv(NOAWT = TRUE) #descomentei para ver + library(rJava) + mx <- maxent(predictors, pres_train) + png(filename = paste0(models.dir, "/", sp, "/maxent_variable_contribution_", + sp, "_", i, ".png")) + plot(mx) + dev.off() + + png(filename = paste0(models.dir, "/", sp, "/maxent_variable_response_", sp, + "_", i, ".png")) + response(mx) + dev.off() + + emx <- evaluate(pres_test, backg_test, mx, predictors) + thresholdmx <- emx@t[which.max(emx@TPR + emx@TNR)] + thmx <- threshold(emx) + mx_TSS <- max(emx@TPR + emx@TNR) - 1 + mx_cont <- predict(mx, predictors, progress = "text") + mx_bin <- mx_cont > thresholdmx + mx_cut <- mx_cont * mx_bin + thmx$AUC <- emx@auc + thmx$TSS <- mx_TSS + thmx$algoritmo <- "maxent" + thmx$partition <- i + row.names(thmx) <- paste(sp, i, "maxent") + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + mx_cont <- cropModel(mx_cont, mask) + mx_bin <- cropModel(mx_bin, mask) + mx_cut <- cropModel(mx_cut, mask) + } + writeRaster(x = mx_cont, filename = paste0(models.dir, "/", sp, "/maxent_cont_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = mx_bin, filename = paste0(models.dir, "/", sp, "/maxent_bin_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = mx_cut, filename = paste0(models.dir, "/", sp, "/maxent_cut_", + sp, "_", i, ".tif"), overwrite = T) + + png(filename = paste0(models.dir, "/", sp, "/maxent", sp, "_", i, "%03d.png")) + plot(mx_cont, main = paste("Maxent raw", "\n", "AUC =", round(emx@auc, 2), "-", + "TSS =", round(mx_TSS, 2))) + plot(mx_bin, main = paste("Maxent P/A", "\n", "AUC =", round(emx@auc, 2), "-", + "TSS =", round(mx_TSS, 2))) + plot(mx_cut, main = paste("Maxent cut", "\n", "AUC =", round(emx@auc, 2), "-", + "TSS =", round(mx_TSS, 2))) + dev.off() + + if (project.model == T) { + for (proj in projections) { + data <- list.files(paste0("./env/", proj), pattern = proj) + data2 <- stack(data) + mx_proj <- predict(data2, mx, progress = "text") + mx_proj_bin <- mx_proj > thresholdmx + mx_proj_cut <- mx_proj_bin * mx_proj + # Normaliza o modelo cut mx_proj_cut <- mx_proj_cut/maxValue(mx_proj_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { source("./fct/cropModel.R") - mx_cont <- cropModel(mx_cont , mask) - mx_bin <- cropModel(mx_bin , mask) - mx_cut <- cropModel(mx_cut , mask) + mx_proj <- cropModel(mx_proj, mask) + mx_proj_bin <- cropModel(mx_proj_bin, mask) + mx_proj_cut <- cropModel(mx_proj_cut, mask) + } + writeRaster(x = mx_proj, filename = paste0(models.dir, "/", sp, "/", + proj, "/maxent_cont_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = mx_proj_bin, filename = paste0(models.dir, "/", sp, "/", + proj, "/maxent_bin_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = mx_proj_cut, filename = paste0(models.dir, "/", sp, "/", + proj, "/maxent_cut_", sp, "_", i, ".tif"), overwrite = T) + rm(data2) } - writeRaster(x=mx_cont,filename=paste0(models.dir,"/",sp,"/maxent_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=mx_bin,filename=paste0(models.dir,"/",sp,"/maxent_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=mx_cut,filename=paste0(models.dir,"/",sp,"/maxent_cut_",sp,"_",i,".tif"),overwrite=T) - - png(filename=paste0(models.dir,"/",sp,"/maxent",sp,"_",i,"%03d.png")) - plot(mx_cont,main=paste("Maxent raw","\n","AUC =", round(emx@auc,2),'-',"TSS =",round(mx_TSS,2))) - plot(mx_bin,main=paste("Maxent P/A","\n","AUC =", round(emx@auc,2),'-',"TSS =",round(mx_TSS,2))) - plot(mx_cut,main=paste("Maxent cut","\n","AUC =", round(emx@auc,2),'-',"TSS =",round(mx_TSS,2))) - dev.off() - - if (project.model == T){ - for (proj in projections){ - data <- list.files(paste0("./env/",proj),pattern=proj) - data2 <- stack(data) - mx_proj <- predict(data2,mx,progress='text') - mx_proj_bin <- mx_proj > thresholdmx - mx_proj_cut <- mx_proj_bin * mx_proj - # Normaliza o modelo cut - #mx_proj_cut <- mx_proj_cut/maxValue(mx_proj_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - mx_proj <- cropModel(mx_proj , mask) - mx_proj_bin <- cropModel(mx_proj_bin , mask) - mx_proj_cut <- cropModel(mx_proj_cut , mask) - } - writeRaster(x=mx_proj,filename=paste0(models.dir,"/",sp,"/",proj,"/maxent_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=mx_proj_bin,filename=paste0(models.dir,"/",sp,"/",proj,"/maxent_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=mx_proj_cut,filename=paste0(models.dir,"/",sp,"/",proj,"/maxent_cut_",sp,"_",i,".tif"),overwrite=T) - rm(data2) - } - } - return(thmx) + } + return(thmx) } do_GLM <- function(sp, - predictors = predictors, - sdmdata_train = sdmdata_train, - envtest_pre = envtest_pre, - envtest_back = envtest_back, - i = i, - envtrain = envtrain, - models.dir = model.dir, - project.model = project.model, - projections = projections, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ){## - cat(paste("GLM",'\n')) - null.model <- glm(sdmdata_train$pa~1,data=envtrain,family="binomial") - full.model <- glm(sdmdata_train$pa~.,data=envtrain,family="binomial") - glm <- step(object = null.model,scope = formula(full.model),direction = "both",trace=F) - eglm <- evaluate(envtest_pre,envtest_back,model=glm,type="response")##### - #eglm <- evaluate(pres_test,backg_test,glm,predictors,type="response") - glm_TSS <- max(eglm@TPR + eglm@TNR)-1 - thresholdglm <- eglm@t[which.max(eglm@TPR + eglm@TNR)] - thglm <- threshold (eglm) - thglm$AUC <- eglm@auc - thglm$TSS <- glm_TSS - thglm$algoritmo <- "glm" - thglm$partition <- i - row.names(thglm) <- paste(sp,i,"glm") - - glm_cont <- predict(predictors,glm,progress='text',type="response") - glm_bin <- glm_cont>thresholdglm - glm_cut <- glm_bin * glm_cont - # Normaliza o modelo cut - #glm_cut <- glm_cut/maxValue(glm_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ + predictors = predictors, + sdmdata_train = sdmdata_train, + envtest_pre = envtest_pre, + envtest_back = envtest_back, + i = i, + envtrain = envtrain, + models.dir = model.dir, + project.model = project.model, + projections = projections, + mask = mask) { + ## + cat(paste("GLM", "\n")) + null.model <- glm(sdmdata_train$pa ~ 1, data = envtrain, family = "binomial") + full.model <- glm(sdmdata_train$pa ~ ., data = envtrain, family = "binomial") + glm <- step(object = null.model, scope = formula(full.model), direction = "both", + trace = F) + eglm <- evaluate(envtest_pre, envtest_back, model = glm, type = "response") ##### + # eglm <- evaluate(pres_test,backg_test,glm,predictors,type='response') + glm_TSS <- max(eglm@TPR + eglm@TNR) - 1 + thresholdglm <- eglm@t[which.max(eglm@TPR + eglm@TNR)] + thglm <- threshold(eglm) + thglm$AUC <- eglm@auc + thglm$TSS <- glm_TSS + thglm$algoritmo <- "glm" + thglm$partition <- i + row.names(thglm) <- paste(sp, i, "glm") + + glm_cont <- predict(predictors, glm, progress = "text", type = "response") + glm_bin <- glm_cont > thresholdglm + glm_cut <- glm_bin * glm_cont + # Normaliza o modelo cut glm_cut <- glm_cut/maxValue(glm_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + glm_cont <- cropModel(glm_cont, mask) + glm_bin <- cropModel(glm_bin, mask) + glm_cut <- cropModel(glm_cut, mask) + } + writeRaster(x = glm_cont, filename = paste0(models.dir, "/", sp, "/glm_cont_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = glm_bin, filename = paste0(models.dir, "/", sp, "/glm_bin_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = glm_cut, filename = paste0(models.dir, "/", sp, "/glm_cut_", + sp, "_", i, ".tif"), overwrite = T) + + png(filename = paste0(models.dir, "/", sp, "/glm", sp, "_", i, "%03d.png")) + plot(glm_cont, main = paste("GLM raw", "\n", "AUC =", round(eglm@auc, 2), "-", + "TSS =", round(glm_TSS, 2))) + plot(glm_bin, main = paste("GLM P/A", "\n", "AUC =", round(eglm@auc, 2), "-", + "TSS =", round(glm_TSS, 2))) + plot(glm_cut, main = paste("GLM cut", "\n", "AUC =", round(eglm@auc, 2), "-", + "TSS =", round(glm_TSS, 2))) + dev.off() + + if (project.model == T) { + for (proj in projections) { + data <- list.files(paste0("./env/", proj), pattern = proj) + data2 <- stack(data) + glm_proj <- predict(data2, glm, progress = "text") + glm_proj_bin <- glm_proj > thresholdglm + glm_proj_cut <- glm_proj_bin * glm_proj + # Normaliza o modelo cut glm_proj_cut <- glm_proj_cut/maxValue(glm_proj_cut) + + if (class(mask) == "SpatialPolygonsDataFrame") { source("./fct/cropModel.R") - glm_cont <- cropModel(glm_cont , mask) - glm_bin <- cropModel(glm_bin , mask) - glm_cut <- cropModel(glm_cut , mask) - } - writeRaster(x=glm_cont,filename=paste0(models.dir,"/",sp,"/glm_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=glm_bin,filename=paste0(models.dir,"/",sp,"/glm_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=glm_cut,filename=paste0(models.dir,"/",sp,"/glm_cut_",sp,"_",i,".tif"),overwrite=T) - - png(filename=paste0(models.dir,"/",sp,"/glm",sp,"_",i,"%03d.png")) - plot(glm_cont,main=paste("GLM raw","\n","AUC =", round(eglm@auc,2),'-',"TSS =",round(glm_TSS,2))) - plot(glm_bin,main=paste("GLM P/A","\n","AUC =", round(eglm@auc,2),'-',"TSS =",round(glm_TSS,2))) - plot(glm_cut,main=paste("GLM cut","\n","AUC =", round(eglm@auc,2),'-',"TSS =",round(glm_TSS,2))) - dev.off() - - if (project.model == T){ - for (proj in projections){ - data <- list.files(paste0("./env/",proj),pattern=proj) - data2 <- stack(data) - glm_proj <- predict(data2,glm,progress='text') - glm_proj_bin <- glm_proj > thresholdglm - glm_proj_cut <- glm_proj_bin * glm_proj - # Normaliza o modelo cut - #glm_proj_cut <- glm_proj_cut/maxValue(glm_proj_cut) - - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - glm_proj <- cropModel(glm_proj , mask) - glm_proj_bin <- cropModel(glm_proj_bin , mask) - glm_proj_cut <- cropModel(glm_proj_cut , mask) - } - writeRaster(x=glm_proj,filename=paste0(models.dir,"/",sp,"/",proj,"/glm_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=glm_proj_bin,filename=paste0(models.dir,"/",sp,"/",proj,"/glm_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=glm_proj_cut,filename=paste0(models.dir,"/",sp,"/",proj,"/glm_cut_",sp,"_",i,".tif"),overwrite=T) - rm(data2) - } + glm_proj <- cropModel(glm_proj, mask) + glm_proj_bin <- cropModel(glm_proj_bin, mask) + glm_proj_cut <- cropModel(glm_proj_cut, mask) + } + writeRaster(x = glm_proj, filename = paste0(models.dir, "/", sp, "/", + proj, "/glm_cont_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = glm_proj_bin, filename = paste0(models.dir, "/", sp, + "/", proj, "/glm_bin_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = glm_proj_cut, filename = paste0(models.dir, "/", sp, + "/", proj, "/glm_cut_", sp, "_", i, ".tif"), overwrite = T) + rm(data2) } - return(thglm) + } + return(thglm) } do_domain <- function(sp, - predictors = predictors, - pres_train = pres_train, #NEW - pres_test = pres_test, #NEW - backg_test = backg_test, #NEW - i = i, #NEW - models.dir = model.dir, - project.model = project.model, - projections = projections, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ){ - cat(paste("Domain",'\n')) - do <- domain (predictors, pres_train) - edo <- evaluate(pres_test,backg_test,do,predictors) - thresholddo <- edo@t[which.max(edo@TPR + edo@TNR)] - thdo <- threshold(edo) - do_TSS <- max(edo@TPR + edo@TNR)-1 - do_cont <- predict(predictors, do, progress='text') - do_bin <- do_cont > thresholddo - do_cut <- do_cont * do_bin - thdo$AUC <- edo@auc - thdo$TSS <- do_TSS - thdo$algoritmo <- "Domain" - thdo$partition <- i - row.names(thdo) <- paste(sp,i,"Domain") - - if (class(mask) == "SpatialPolygonsDataFrame"){ + predictors = predictors, + pres_train = pres_train, + pres_test = pres_test, + backg_test = backg_test, + i = i, + models.dir = model.dir, + project.model = project.model, + projections = projections, + mask = mask) { + cat(paste("Domain", "\n")) + do <- domain(predictors, pres_train) + edo <- evaluate(pres_test, backg_test, do, predictors) + thresholddo <- edo@t[which.max(edo@TPR + edo@TNR)] + thdo <- threshold(edo) + do_TSS <- max(edo@TPR + edo@TNR) - 1 + do_cont <- predict(predictors, do, progress = "text") + do_bin <- do_cont > thresholddo + do_cut <- do_cont * do_bin + thdo$AUC <- edo@auc + thdo$TSS <- do_TSS + thdo$algoritmo <- "Domain" + thdo$partition <- i + row.names(thdo) <- paste(sp, i, "Domain") + + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + do_cont <- cropModel(do_cont, mask) + do_bin <- cropModel(do_bin, mask) + do_cut <- cropModel(do_cut, mask) + } + writeRaster(x = do_cont, filename = paste0(models.dir, "/", sp, "/Domain_cont_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = do_bin, filename = paste0(models.dir, "/", sp, "/Domain_bin_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = do_cut, filename = paste0(models.dir, "/", sp, "/Domain_cut_", + sp, "_", i, ".tif"), overwrite = T) + + png(filename = paste0(models.dir, "/", sp, "/Domain", sp, "_", i, "%03d.png")) + plot(do_cont, main = paste("Domain raw", "\n", "AUC =", round(edo@auc, 2), "-", + "TSS =", round(do_TSS, 2))) + plot(do_bin, main = paste("Domain P/A", "\n", "AUC =", round(edo@auc, 2), "-", + "TSS =", round(do_TSS, 2))) + plot(do_cut, main = paste("Domain cut", "\n", "AUC =", round(edo@auc, 2), "-", + "TSS =", round(do_TSS, 2))) + dev.off() + + + if (project.model == T) { + for (proj in projections) { + data <- list.files(paste0("./env/", proj), pattern = proj) + data2 <- stack(data) + do_proj <- predict(data2, do, progress = "text") + do_proj_bin <- do_proj > thresholddo + do_proj_cut <- do_proj_bin * do_proj + # Normaliza o modelo cut do_proj_cut <- do_proj_cut/maxValue(do_proj_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { source("./fct/cropModel.R") - do_cont <- cropModel(do_cont , mask) - do_bin <- cropModel(do_bin , mask) - do_cut <- cropModel(do_cut , mask) - } - writeRaster(x=do_cont,filename=paste0(models.dir,"/",sp,"/Domain_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=do_bin,filename=paste0(models.dir,"/",sp,"/Domain_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=do_cut,filename=paste0(models.dir,"/",sp,"/Domain_cut_",sp,"_",i,".tif"),overwrite=T) - - png(filename=paste0(models.dir,"/",sp,"/Domain",sp,"_",i,"%03d.png")) - plot(do_cont,main=paste("Domain raw","\n","AUC =", round(edo@auc,2),'-',"TSS =",round(do_TSS,2))) - plot(do_bin,main=paste("Domain P/A","\n","AUC =", round(edo@auc,2),'-',"TSS =",round(do_TSS,2))) - plot(do_cut,main=paste("Domain cut","\n","AUC =", round(edo@auc,2),'-',"TSS =",round(do_TSS,2))) - dev.off() - - - if (project.model == T){ - for (proj in projections){ - data <- list.files(paste0("./env/",proj),pattern=proj) - data2 <- stack(data) - do_proj <- predict(data2,do,progress='text') - do_proj_bin <- do_proj > thresholddo - do_proj_cut <- do_proj_bin * do_proj - # Normaliza o modelo cut - #do_proj_cut <- do_proj_cut/maxValue(do_proj_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - do_proj <- cropModel(do_proj , mask) - do_proj_bin <- cropModel(do_proj_bin , mask) - do_proj_cut <- cropModel(do_proj_cut , mask) - } - writeRaster(x=do_proj,filename=paste0(models.dir,"/",sp,"/",proj,"/Domain_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=do_proj_bin,filename=paste0(models.dir,"/",sp,"/",proj,"/Domain_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=do_proj_cut,filename=paste0(models.dir,"/",sp,"/",proj,"/Domain_cut_",sp,"_",i,".tif"),overwrite=T) - rm(data2) - } + do_proj <- cropModel(do_proj, mask) + do_proj_bin <- cropModel(do_proj_bin, mask) + do_proj_cut <- cropModel(do_proj_cut, mask) + } + writeRaster(x = do_proj, filename = paste0(models.dir, "/", sp, "/", + proj, "/Domain_cont_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = do_proj_bin, filename = paste0(models.dir, "/", sp, "/", + proj, "/Domain_bin_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = do_proj_cut, filename = paste0(models.dir, "/", sp, "/", + proj, "/Domain_cut_", sp, "_", i, ".tif"), overwrite = T) + rm(data2) } - return(thdo) + } + return(thdo) } do_mahal <- function(sp, - predictors = predictors, - pres_train = pres_train, #NEW - pres_test = pres_test, #NEW - backg_test = backg_test, #NEW - i = i, #NEW - models.dir = model.dir, - project.model = project.model, - projections = projections, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ){ - cat(paste("Mahalanobis distance",'\n')) - ma <- mahal (predictors, pres_train) - if (exists("ma")){ - ema <- evaluate(pres_test,backg_test,ma,predictors) - thresholdma <- ema@t[which.max(ema@TPR + ema@TNR)] - thma <- threshold(ema) - ma_TSS <- max(ema@TPR + ema@TNR)-1 - ma_cont <- predict(ma, predictors,progress='text') - ma_cont[ma_cont < threshold(ema,'no_omission')] <- threshold(ema,'no_omission') - ma_bin <- ma_cont > thresholdma - ma_cut <- ma_cont - ma_cut[ma_cut < thresholdma] <- thresholdma - if(minValue(ma_cut)<0) { - ma_cut<-(ma_cut-minValue(ma_cut))/maxValue(ma_cut-minValue(ma_cut))} - - thma$AUC <- ema@auc - thma$TSS <- ma_TSS - thma$algoritmo <- "Mahal" - thma$partition <- i - row.names(thma) <- paste(sp,i,"Mahal") - - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - ma_cont <- cropModel(ma_cont , mask) - ma_bin <- cropModel(ma_bin , mask) - ma_cut <- cropModel(ma_cut , mask) - } - writeRaster(x=ma_cont,filename=paste0(models.dir,"/",sp,"/Mahal_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=ma_bin,filename=paste0(models.dir,"/",sp,"/Mahal_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=ma_cut,filename=paste0(models.dir,"/",sp,"/Mahal_cut_",sp,"_",i,".tif"),overwrite=T) - - png(filename=paste0(models.dir,"/",sp,"/Mahal",sp,"_",i,"%03d.png")) - plot(ma_cont,main=paste("Mahal raw","\n","AUC =", round(ema@auc,2),'-',"TSS =",round(ma_TSS,2))) - plot(ma_bin,main=paste("Mahal P/A","\n","AUC =", round(ema@auc,2),'-',"TSS =",round(ma_TSS,2))) - plot(ma_cut,main=paste("Mahal cut","\n","AUC =", round(ema@auc,2),'-',"TSS =",round(ma_TSS,2))) - dev.off() - - if (project.model == T){ - for (proj in projections){ - data <- list.files(paste0("./env/",proj),pattern=proj) - data2 <- stack(data) - ma_proj <- predict(data2,ma,progress='text') - ma_proj_bin <- ma_proj > thresholdma - ma_proj_cut <- ma_proj_bin * ma_proj - # Normaliza o modelo cut - #ma_proj_cut <- ma_proj_cut/maxValue(ma_proj_cut) - - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - ma_proj <- cropModel(ma_proj , mask) - ma_proj_bin <- cropModel(ma_proj_bin , mask) - ma_proj_cut <- cropModel(ma_proj_cut , mask) - } - writeRaster(x=ma_proj,filename=paste0(models.dir,"/",sp,"/",proj,"/mahal_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=ma_proj_bin,filename=paste0(models.dir,"/",sp,"/",proj,"/mahal_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=ma_proj_cut,filename=paste0(models.dir,"/",sp,"/",proj,"/mahal_cut_",sp,"_",i,".tif"),overwrite=T) - rm(data2) - } + predictors = predictors, + pres_train = pres_train, + pres_test = pres_test, + backg_test = backg_test, + i = i, + models.dir = model.dir, + project.model = project.model, + projections = projections, + mask = mask) { + cat(paste("Mahalanobis distance", "\n")) + ma <- mahal(predictors, pres_train) + if (exists("ma")) { + ema <- evaluate(pres_test, backg_test, ma, predictors) + thresholdma <- ema@t[which.max(ema@TPR + ema@TNR)] + thma <- threshold(ema) + ma_TSS <- max(ema@TPR + ema@TNR) - 1 + ma_cont <- predict(ma, predictors, progress = "text") + ma_cont[ma_cont < threshold(ema, "no_omission")] <- threshold(ema, "no_omission") + ma_bin <- ma_cont > thresholdma + ma_cut <- ma_cont + ma_cut[ma_cut < thresholdma] <- thresholdma + if (minValue(ma_cut) < 0) { + ma_cut <- (ma_cut - minValue(ma_cut))/maxValue(ma_cut - minValue(ma_cut)) + } + + thma$AUC <- ema@auc + thma$TSS <- ma_TSS + thma$algoritmo <- "Mahal" + thma$partition <- i + row.names(thma) <- paste(sp, i, "Mahal") + + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + ma_cont <- cropModel(ma_cont, mask) + ma_bin <- cropModel(ma_bin, mask) + ma_cut <- cropModel(ma_cut, mask) + } + writeRaster(x = ma_cont, filename = paste0(models.dir, "/", sp, "/Mahal_cont_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = ma_bin, filename = paste0(models.dir, "/", sp, "/Mahal_bin_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = ma_cut, filename = paste0(models.dir, "/", sp, "/Mahal_cut_", + sp, "_", i, ".tif"), overwrite = T) + + png(filename = paste0(models.dir, "/", sp, "/Mahal", sp, "_", i, "%03d.png")) + plot(ma_cont, main = paste("Mahal raw", "\n", "AUC =", round(ema@auc, 2), + "-", "TSS =", round(ma_TSS, 2))) + plot(ma_bin, main = paste("Mahal P/A", "\n", "AUC =", round(ema@auc, 2), + "-", "TSS =", round(ma_TSS, 2))) + plot(ma_cut, main = paste("Mahal cut", "\n", "AUC =", round(ema@auc, 2), + "-", "TSS =", round(ma_TSS, 2))) + dev.off() + + if (project.model == T) { + for (proj in projections) { + data <- list.files(paste0("./env/", proj), pattern = proj) + data2 <- stack(data) + ma_proj <- predict(data2, ma, progress = "text") + ma_proj_bin <- ma_proj > thresholdma + ma_proj_cut <- ma_proj_bin * ma_proj + # Normaliza o modelo cut ma_proj_cut <- ma_proj_cut/maxValue(ma_proj_cut) + + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + ma_proj <- cropModel(ma_proj, mask) + ma_proj_bin <- cropModel(ma_proj_bin, mask) + ma_proj_cut <- cropModel(ma_proj_cut, mask) } + writeRaster(x = ma_proj, filename = paste0(models.dir, "/", sp, "/", + proj, "/mahal_cont_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = ma_proj_bin, filename = paste0(models.dir, "/", sp, + "/", proj, "/mahal_bin_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = ma_proj_cut, filename = paste0(models.dir, "/", sp, + "/", proj, "/mahal_cut_", sp, "_", i, ".tif"), overwrite = T) + rm(data2) + } } - else cat("Mahalanobis distance did not run") - return(thma) + } else cat("Mahalanobis distance did not run") + return(thma) } do_SVM2 <- function(sp, - predictors = predictors, - sdmdata_train = sdmdata_train, - envtest_pre = envtest_pre, - envtest_back = envtest_back, - i = i, - envtrain = envtrain, - models.dir = model.dir, - project.model = project.model, - projections = projections, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ){ - cat(paste("SVM2",'\n')) - library(e1071) - svm2 <- best.tune("svm",envtrain,sdmdata_train$pa,data=envtrain)##svm deve ser com a variável resposta binária ou contínua, eu acho que binária - esvm2 <- evaluate(envtest_pre,envtest_back,svm2) - #esvm <- evaluate(pres_test,backg_test,model = svm,x = predictors) - svm2_TSS <- max(esvm2@TPR + esvm2@TNR)-1 - thresholdsvm2 <- esvm2@t[which.max(esvm2@TPR + esvm2@TNR)] - thsvm2 <- threshold (esvm2) - thsvm2$AUC <- esvm2@auc - thsvm2$TSS <- svm2_TSS - thsvm2$algoritmo <- "svm2" - thsvm2$partition <- i - row.names(thsvm2) <- paste(sp,i,"svm2") - svm2_cont <- predict(predictors,svm2,progress='text') - svm2_bin <- svm2_cont>thresholdsvm2 - svm2_cut <- svm2_bin * svm2_cont - - #TRANSFORMA 0 A 1 - svm2_cont <- svm2_cont/maxValue(svm2_cont) - svm2_cut <- svm2_cut/maxValue(svm2_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ + predictors = predictors, + sdmdata_train = sdmdata_train, + envtest_pre = envtest_pre, + envtest_back = envtest_back, + i = i, + envtrain = envtrain, + models.dir = model.dir, + project.model = project.model, + projections = projections, + mask = mask) { + cat(paste("SVM2", "\n")) + library(e1071) + svm2 <- best.tune("svm", envtrain, sdmdata_train$pa, data = envtrain) ##svm deve ser com a variável resposta binária ou contínua, eu acho que binária + esvm2 <- evaluate(envtest_pre, envtest_back, svm2) + # esvm <- evaluate(pres_test,backg_test,model = svm,x = predictors) + svm2_TSS <- max(esvm2@TPR + esvm2@TNR) - 1 + thresholdsvm2 <- esvm2@t[which.max(esvm2@TPR + esvm2@TNR)] + thsvm2 <- threshold(esvm2) + thsvm2$AUC <- esvm2@auc + thsvm2$TSS <- svm2_TSS + thsvm2$algoritmo <- "svm2" + thsvm2$partition <- i + row.names(thsvm2) <- paste(sp, i, "svm2") + svm2_cont <- predict(predictors, svm2, progress = "text") + svm2_bin <- svm2_cont > thresholdsvm2 + svm2_cut <- svm2_bin * svm2_cont + + # TRANSFORMA 0 A 1 + svm2_cont <- svm2_cont/maxValue(svm2_cont) + svm2_cut <- svm2_cut/maxValue(svm2_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { + source("./fct/cropModel.R") + svm2_cont <- cropModel(svm2_cont, mask) + svm2_bin <- cropModel(svm2_bin, mask) + svm2_cut <- cropModel(svm2_cut, mask) + } + writeRaster(x = svm2_cont, filename = paste0(models.dir, "/", sp, "/svm2_cont_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = svm2_bin, filename = paste0(models.dir, "/", sp, "/svm2_bin_", + sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = svm2_cut, filename = paste0(models.dir, "/", sp, "/svm2_cut_", + sp, "_", i, ".tif"), overwrite = T) + + png(filename = paste0(models.dir, "/", sp, "/svm2", sp, "_", i, "%03d.png")) + plot(svm2_cont, main = paste("SVM2 raw", "\n", "AUC =", round(esvm2@auc, 2), + "-", "TSS =", round(svm2_TSS, 2))) + plot(svm2_bin, main = paste("SVM2 P/A", "\n", "AUC =", round(esvm2@auc, 2), "-", + "TSS =", round(svm2_TSS, 2))) + plot(svm2_cut, main = paste("SVM2 cut", "\n", "AUC =", round(esvm2@auc, 2), "-", + "TSS =", round(svm2_TSS, 2))) + dev.off() + + if (project.model == T) { + for (proj in projections) { + data <- list.files(paste0("./env/", proj), pattern = proj) + data2 <- stack(data) + svm2_proj <- predict(data2, svm2, progress = "text") + svm2_proj_bin <- svm2_proj > thresholdsvm2 + svm2_proj_cut <- svm2_proj_bin * svm2_proj + # Normaliza o modelo cut svm2_proj_cut <- svm2_proj_cut/maxValue(svm2_proj_cut) + if (class(mask) == "SpatialPolygonsDataFrame") { source("./fct/cropModel.R") - svm2_cont <- cropModel(svm2_cont , mask) - svm2_bin <- cropModel(svm2_bin , mask) - svm2_cut <- cropModel(svm2_cut , mask) + svm2_proj <- cropModel(svm2_proj, mask) + svm2_proj_bin <- cropModel(svm2_proj_bin, mask) + svm2_proj_cut <- cropModel(svm2_proj_cut, mask) + } + writeRaster(x = svm2_proj, filename = paste0(models.dir, "/", sp, "/", + proj, "/svm2_cont_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = svm2_proj_bin, filename = paste0(models.dir, "/", sp, + "/", proj, "/svm2_bin_", sp, "_", i, ".tif"), overwrite = T) + writeRaster(x = svm2_proj_cut, filename = paste0(models.dir, "/", sp, + "/", proj, "/svm2_cut_", sp, "_", i, ".tif"), overwrite = T) + rm(data2) } - writeRaster(x=svm2_cont,filename=paste0(models.dir,"/",sp,"/svm2_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=svm2_bin,filename=paste0(models.dir,"/",sp,"/svm2_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=svm2_cut,filename=paste0(models.dir,"/",sp,"/svm2_cut_",sp,"_",i,".tif"),overwrite=T) - - png(filename=paste0(models.dir,"/",sp,"/svm2",sp,"_",i,"%03d.png")) - plot(svm2_cont,main=paste("SVM2 raw","\n","AUC =", round(esvm2@auc,2),'-',"TSS =",round(svm2_TSS,2))) - plot(svm2_bin,main=paste("SVM2 P/A","\n","AUC =", round(esvm2@auc,2),'-',"TSS =",round(svm2_TSS,2))) - plot(svm2_cut,main=paste("SVM2 cut","\n","AUC =", round(esvm2@auc,2),'-',"TSS =",round(svm2_TSS,2))) - dev.off() - - if (project.model == T){ - for (proj in projections){ - data <- list.files(paste0("./env/",proj),pattern=proj) - data2 <- stack(data) - svm2_proj <- predict(data2,svm2,progress='text') - svm2_proj_bin <- svm2_proj > thresholdsvm2 - svm2_proj_cut <- svm2_proj_bin * svm2_proj - # Normaliza o modelo cut - #svm2_proj_cut <- svm2_proj_cut/maxValue(svm2_proj_cut) - if (class(mask) == "SpatialPolygonsDataFrame"){ - source("./fct/cropModel.R") - svm2_proj <- cropModel(svm2_proj , mask) - svm2_proj_bin <- cropModel(svm2_proj_bin , mask) - svm2_proj_cut <- cropModel(svm2_proj_cut , mask) - } - writeRaster(x=svm2_proj,filename=paste0(models.dir,"/",sp,"/",proj,"/svm2_cont_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=svm2_proj_bin,filename=paste0(models.dir,"/",sp,"/",proj,"/svm2_bin_",sp,"_",i,".tif"),overwrite=T) - writeRaster(x=svm2_proj_cut,filename=paste0(models.dir,"/",sp,"/",proj,"/svm2_cut_",sp,"_",i,".tif"),overwrite=T) - rm(data2) - } - } - return(thsvm2) + } + return(thsvm2) } -dismo.mod <- function(sp, - occs = spp.filt,#complete occurrence table - predictors = predictors, - buffer = TRUE, - buffer.type = "max",#"mean" - maxent = F, - Bioclim = T, - Domain = F, - Mahal = F, - GLM = F, - RF = F, - SVM = F, - SVM2 = F, - part = 3, - seed = NULL,#for reproducibility purposes - models.dir = "./models", - project.model = F, - projections = NULL, - projdata = NULL,#um vector con nombres - #stack_gcms = "future_vars", # Lista dos stacks de cada GCM. Ex: stack1 <- stack(variaveis_HADGEM); stack2<-stack(variaveis_CANESM); stack_gcms<-c(stack1,stack2) - mask = mask,# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - n.back = 500){ +# dismo.mod() possuia o seguinte parâmetro comentado: stack_gcms = 'future_vars', +# # Lista dos stacks de cada GCM. Ex: stack1 <- stack(variaveis_HADGEM); +# stack2<-stack(variaveis_CANESM); stack_gcms<-c(stack1,stack2) - if (file.exists(paste0(models.dir)) == FALSE) dir.create(paste0(models.dir)) - if (file.exists(paste0(models.dir,"/",sp)) == FALSE) dir.create(paste0(models.dir,"/",sp)) +dismo.mod <- function(sp, + occs = spp.filt, + predictors = predictors, + buffer = TRUE, + buffer.type = "max", + maxent = F, + Bioclim = T, + Domain = F, + Mahal = F, + GLM = F, + RF = F, + SVM = F, + SVM2 = F, + part = 3, + seed = NULL, + models.dir = "./models", + project.model = F, + projections = NULL, + projdata = NULL, + mask = mask, + n.back = 500) { + + if (file.exists(paste0(models.dir)) == FALSE) + dir.create(paste0(models.dir)) + if (file.exists(paste0(models.dir, "/", sp)) == FALSE) + dir.create(paste0(models.dir, "/", sp)) if (project.model == T) { - for (proj in projections){ - if (file.exists(paste0(models.dir,"/",sp,"/",proj)) == FALSE) dir.create(paste0(models.dir,"/",sp,"/",proj)) + for (proj in projections) { + if (file.exists(paste0(models.dir, "/", sp, "/", proj)) == FALSE) + dir.create(paste0(models.dir, "/", sp, "/", proj)) } } - - #Teste se pacotes estao instalados. Caso negativo, instala os pacotes e dependencias - pacotes <- c('dismo','XML','raster','rgdal','maps','rgeos') - # if (!pacotes %in% installed.packages()) install.packages(pacotes, dependencies = TRUE) - #Loading packages + + # Teste se pacotes estao instalados. Caso negativo, instala os pacotes e + # dependencias + pacotes <- c("dismo", "XML", "raster", "rgdal", "maps", "rgeos") + # if (!pacotes %in% installed.packages()) install.packages(pacotes, dependencies + # = TRUE) Loading packages lapply(pacotes, require, character.only = TRUE) - - + + print(date()) - - cat(paste("Modeling",sp,"...",'\n')) - #extrai as coordenadas de cada especie - coord <- occs[occs$sp == sp,c('lon','lat')] + + cat(paste("Modeling", sp, "...", "\n")) + # extrai as coordenadas de cada especie + coord <- occs[occs$sp == sp, c("lon", "lat")] n <- nrow(coord) - n.back <- n * 2 #descomentado - #tabela de valores + n.back <- n * 2 #descomentado + # tabela de valores presvals <- raster::extract(predictors, coord) - - if (buffer == FALSE) { - set.seed(seed+2) - backgr <- randomPoints(predictors, n.back) - } else { - source("./fct/createBuffer.R") -backgr <- createBuffer(coord_ = coord, n.back_ = n.back, buffer.type_ = buffer.type, occs_ = occs, sp_ = sp, seed_ = seed) + + if (buffer == FALSE) { + set.seed(seed + 2) + backgr <- randomPoints(predictors, n.back) + } else { + source("./fct/createBuffer.R") + backgr <- createBuffer(coord_ = coord, n.back_ = n.back, buffer.type_ = buffer.type, + occs_ = occs, sp_ = sp, seed_ = seed) } - - colnames(backgr) <- c('lon', 'lat') - - #Extraindo dados ambientais dos bckgr + + colnames(backgr) <- c("lon", "lat") + + # Extraindo dados ambientais dos bckgr backvals <- raster::extract(predictors, backgr) pa <- c(rep(1, nrow(presvals)), rep(0, nrow(backvals))) - - #Data partition - if (n<11) part <- n else part <- part - set.seed(seed)#reproducibility - group <- kfold(coord,part) - set.seed(seed+1) - bg.grp <- kfold(backgr,part) - group.all <- c(group,bg.grp) - - pres <- cbind(coord,presvals) - back <- cbind(backgr,backvals) - rbind_1 <- rbind(pres,back) - sdmdata <- data.frame(cbind(group.all,pa,rbind_1)) - rm(rbind_1) - rm(pres) - rm(back) - gc() - write.table(sdmdata,file = paste0(models.dir,"/",sp,"/sdmdata.txt")) - - - #####Hace los modelos - for (i in unique(group)){ - cat(paste(sp,"partition number",i,'\n')) + + # Data partition + if (n < 11) + part <- n else part <- part + set.seed(seed) #reproducibility + group <- kfold(coord, part) + set.seed(seed + 1) + bg.grp <- kfold(backgr, part) + group.all <- c(group, bg.grp) + + pres <- cbind(coord, presvals) + back <- cbind(backgr, backvals) + rbind_1 <- rbind(pres, back) + sdmdata <- data.frame(cbind(group.all, pa, rbind_1)) + rm(rbind_1) + rm(pres) + rm(back) + gc() + write.table(sdmdata, file = paste0(models.dir, "/", sp, "/sdmdata.txt")) + + + ##### Hace los modelos + for (i in unique(group)) { + cat(paste(sp, "partition number", i, "\n")) pres_train <- coord[group != i, ] - if(n == 1) pres_train <- coord[group == i,] - pres_test <- coord[group == i, ] - - backg_train <- backgr[bg.grp != i,]#not used? - backg_test <- backgr[bg.grp == i,]#new - - sdmdata_train <- subset(sdmdata,group!=i)#new - sdmdata_test <- subset(sdmdata,group == i)#new - - envtrain <- subset(sdmdata_train,select= c(-group,-lon,-lat))#new - envtest <- subset(sdmdata_test,select=c(-group,-lon,-lat)) - envtest_pre <- subset(sdmdata_test,pa == 1,select= c(-group,-lon,-lat,-pa))#new - envtest_back <- subset(sdmdata_test,pa == 0,select= c(-group,-lon,-lat,-pa))#new - + if (n == 1) + pres_train <- coord[group == i, ] + pres_test <- coord[group == i, ] + + backg_train <- backgr[bg.grp != i, ] #not used? + backg_test <- backgr[bg.grp == i, ] #new + + sdmdata_train <- subset(sdmdata, group != i) #new + sdmdata_test <- subset(sdmdata, group == i) #new + + envtrain <- subset(sdmdata_train, select = c(-group, -lon, -lat)) #new + envtest <- subset(sdmdata_test, select = c(-group, -lon, -lat)) + envtest_pre <- subset(sdmdata_test, pa == 1, select = c(-group, -lon, -lat, + -pa)) #new + envtest_back <- subset(sdmdata_test, pa == 0, select = c(-group, -lon, -lat, + -pa)) #new + ##### Creates a .png plot of the initial dataset - cat(paste("Plotting the dataset...",'\n')) + cat(paste("Plotting the dataset...", "\n")) extent <- extent(predictors) - png(filename=paste0(models.dir,"/",sp,"/",i,sp,"dataset.png")) - par(mfrow=c(1,1),mar=c(5,4,3,0)) - plot(predictors[[1]]!=0,col="grey95",main=paste(sp,"part.",i),legend=F) - map('world',c('',"South America"),xlim=c(extent@xmin,extent@xmax),ylim=c(extent@ymin,extent@ymax),add=T) - points(pres_train,pch=21,bg="red") - points(pres_test,pch=21,bg="blue") - points(backg_train,pch=3,col="red",cex=0.3) - points(backg_test,pch=3,col="blue",cex=0.3) - legend("topright",pch=c(21,21,3,3),pt.bg=c("red","blue"),legend=c("PresTrain","PresTest","BackTrain","BackTest"),col=c("black","black","red","blue")) + png(filename = paste0(models.dir, "/", sp, "/", i, sp, "dataset.png")) + par(mfrow = c(1, 1), mar = c(5, 4, 3, 0)) + plot(predictors[[1]] != 0, col = "grey95", main = paste(sp, "part.", i), + legend = F) + map("world", c("", "South America"), xlim = c(extent@xmin, extent@xmax), + ylim = c(extent@ymin, extent@ymax), add = T) + points(pres_train, pch = 21, bg = "red") + points(pres_test, pch = 21, bg = "blue") + points(backg_train, pch = 3, col = "red", cex = 0.3) + points(backg_test, pch = 3, col = "blue", cex = 0.3) + legend("topright", pch = c(21, 21, 3, 3), pt.bg = c("red", "blue"), legend = c("PresTrain", + "PresTest", "BackTrain", "BackTest"), col = c("black", "black", "red", + "blue")) dev.off() - - #SSB <- ssb(pres_test,backg_train,pres_train) - #cat("Spatial sorting bias (approaches 0 with bias)") - #print(SSB[,1]/SSB[,2]) - #faz os modelos - - cat(paste("Modeling",sp,"Partition",i,'\n')) - eval_df <- data.frame("kappa"=1,"spec_sens"=1,"no_omission"=1,"prevalence"=1,"equal_sens_spec"=1,"sensitivity"=1,"AUC"=1,"TSS"=1,"algoritmo"="foo","partition"=1) - - if (Bioclim == T){ - thbc <- do_bioclim(sp = sp, - predictors = predictors, - pres_train = pres_train, - pres_test = pres_test, - backg_test = backg_test, - i = i, - models.dir = models.dir, - project.model = F, - projections = NULL, - projdata = projdata, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ) - eval_df <- rbind(eval_df, thbc) + + # SSB <- ssb(pres_test,backg_train,pres_train) cat('Spatial sorting bias + # (approaches 0 with bias)') print(SSB[,1]/SSB[,2]) faz os modelos + + cat(paste("Modeling", sp, "Partition", i, "\n")) + eval_df <- data.frame(kappa = 1, spec_sens = 1, no_omission = 1, prevalence = 1, + equal_sens_spec = 1, sensitivity = 1, AUC = 1, TSS = 1, algoritmo = "foo", + partition = 1) + + if (Bioclim == T) { + thbc <- do_bioclim(sp = sp, predictors = predictors, pres_train = pres_train, + pres_test = pres_test, backg_test = backg_test, i = i, models.dir = models.dir, + project.model = F, projections = NULL, projdata = projdata, mask = mask) + eval_df <- rbind(eval_df, thbc) } - - if (Domain == T){ - thdo <- do_domain(sp = sp, - predictors = predictors, - pres_train = pres_train, - pres_test = pres_test, - backg_test = backg_test, - i = i, - models.dir = models.dir, - project.model = F, - projections = NULL, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ) - - eval_df <- rbind(eval_df, thdo) + + if (Domain == T) { + thdo <- do_domain(sp = sp, predictors = predictors, pres_train = pres_train, + pres_test = pres_test, backg_test = backg_test, i = i, models.dir = models.dir, + project.model = F, projections = NULL, mask = mask) + + eval_df <- rbind(eval_df, thdo) } - - if (maxent == T){ - thmx <- do_maxent(sp = sp, - predictors = predictors, - pres_train = pres_train, - pres_test = pres_test, - backg_test = backg_test, - i = i, - models.dir = models.dir, - project.model = F, - projections = NULL, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ) - eval_df <- rbind(eval_df, thmx) + + if (maxent == T) { + thmx <- do_maxent(sp = sp, predictors = predictors, pres_train = pres_train, + pres_test = pres_test, backg_test = backg_test, i = i, models.dir = models.dir, + project.model = F, projections = NULL, mask = mask) + eval_df <- rbind(eval_df, thmx) } - - if (Mahal == T){ - thma <- do_mahal(sp = sp, - predictors = predictors, - pres_train = pres_train, - pres_test = pres_test, - backg_test = backg_test, - i = i, - models.dir = models.dir, - project.model = F, - projections = NULL, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ) - eval_df <- rbind(eval_df, thma) - } - - if (GLM == T){ - thglm <- do_GLM (sp, - predictors = predictors, - sdmdata_train = sdmdata_train, - envtest_pre = envtest_pre, - envtest_back = envtest_back, - i = i, - envtrain = envtrain, - models.dir = models.dir, - project.model = F, - projections = NULL, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ) - eval_df <- rbind(eval_df, thglm) + + if (Mahal == T) { + thma <- do_mahal(sp = sp, predictors = predictors, pres_train = pres_train, + pres_test = pres_test, backg_test = backg_test, i = i, models.dir = models.dir, + project.model = F, projections = NULL, mask = mask) + eval_df <- rbind(eval_df, thma) } - - if (RF == T){ - thrf <- do_randomForest(sp, - predictors = predictors, - sdmdata_train = sdmdata_train, - envtest_pre = envtest_pre, - envtest_back = envtest_back, - i = i, - envtrain = envtrain, - models.dir = models.dir, - project.model = F, - projections = NULL, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ) - eval_df <- rbind(eval_df, thrf) + + if (GLM == T) { + thglm <- do_GLM(sp, predictors = predictors, sdmdata_train = sdmdata_train, + envtest_pre = envtest_pre, envtest_back = envtest_back, i = i, envtrain = envtrain, + models.dir = models.dir, project.model = F, projections = NULL, mask = mask) + eval_df <- rbind(eval_df, thglm) } - - if (SVM == T){ - thsvm <- do_SVM (sp, - predictors = predictors, - sdmdata_train = sdmdata_train, - envtest_pre = envtest_pre, - envtest_back = envtest_back, - i = i, - envtrain = envtrain, - part = part, - models.dir = models.dir, - project.model = F, - projections = NULL, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ) - eval_df <- rbind(eval_df, thsvm) + + if (RF == T) { + thrf <- do_randomForest(sp, predictors = predictors, sdmdata_train = sdmdata_train, + envtest_pre = envtest_pre, envtest_back = envtest_back, i = i, envtrain = envtrain, + models.dir = models.dir, project.model = F, projections = NULL, mask = mask) + eval_df <- rbind(eval_df, thrf) } - - if (SVM2 == T){ - thsvm2 <- do_SVM2(sp, - predictors = predictors, - sdmdata_train = sdmdata_train, - envtest_pre = envtest_pre, - envtest_back = envtest_back, - i = i, - envtrain = envtrain, - models.dir = models.dir, - project.model = F, - projections = NULL, - mask = mask# a SpatialPolygonsDataFrame layer to mask and crop the predicted model - ) - eval_df <- rbind(eval_df, thsvm2) + + if (SVM == T) { + thsvm <- do_SVM(sp, predictors = predictors, sdmdata_train = sdmdata_train, + envtest_pre = envtest_pre, envtest_back = envtest_back, i = i, envtrain = envtrain, + part = part, models.dir = models.dir, project.model = F, projections = NULL, + mask = mask) + eval_df <- rbind(eval_df, thsvm) } - - cat(paste("Saving the evaluation file...",sp,i,'\n')) - write.table(eval_df[-1,],file = paste0(models.dir,"/",sp,"/evaluate",sp,"_",i,".txt")) + + if (SVM2 == T) { + thsvm2 <- do_SVM2(sp, predictors = predictors, sdmdata_train = sdmdata_train, + envtest_pre = envtest_pre, envtest_back = envtest_back, i = i, envtrain = envtrain, + models.dir = models.dir, project.model = F, projections = NULL, mask = mask) + eval_df <- rbind(eval_df, thsvm2) } - - #cat(paste("Saving the evaluation file...",sp,i,'\n')) - #write.table(eval,file = paste0(models.dir,"/",sp,"/evaluate",sp,"_",i,".txt")) - cat(as.character(sp), "DONE on", date(), '\n') - print(date()) + + cat(paste("Saving the evaluation file...", sp, i, "\n")) + write.table(eval_df[-1, ], file = paste0(models.dir, "/", sp, "/evaluate", + sp, "_", i, ".txt")) } + + # cat(paste('Saving the evaluation file...',sp,i,'\n')) write.table(eval,file = + # paste0(models.dir,'/',sp,'/evaluate',sp,'_',i,'.txt')) + cat(as.character(sp), "DONE on", date(), "\n") + print(date()) +} diff --git a/ENM/fct/read.eval1.R b/ENM/fct/read.eval1.R index c80c773..87891d0 100644 --- a/ENM/fct/read.eval1.R +++ b/ENM/fct/read.eval1.R @@ -1,116 +1,128 @@ #####----- -#Function to organize the evaluate output from SDM to be used with xonsolidate.data(). -read.eval<-function(sp, - input.folder1="./FLORA_buffermax", - input.folder2="./FLORA_sem_buffer", - output.folder="./TSS_Evaluate"){ - -library("data.table") +# Function to organize the evaluate output from SDM to be used with +# xonsolidate.data(). +read.eval <- function(sp, input.folder1 = "./FLORA_buffermax", input.folder2 = "./FLORA_sem_buffer", + output.folder = "./TSS_Evaluate") { + + library("data.table") - evall.buffer <- list.files(path = paste0(input.folder1,"/",sp), pattern=paste0("evaluate",sp),full.names = T) - evall <- list.files(path = paste0(input.folder2,"/",sp), pattern=paste0("evaluate",sp),full.names = T) + evall.buffer <- list.files(path = paste0(input.folder1, "/", sp), pattern = paste0("evaluate", + sp), full.names = T) + evall <- list.files(path = paste0(input.folder2, "/", sp), pattern = paste0("evaluate", + sp), full.names = T) # Resultados com Buffer ---- - lista.buffer<-list() - lista.buffer.final<-list() + lista.buffer <- list() + lista.buffer.final <- list() for (i in 1:length(evall.buffer)) { - lista.buffer[[i]]<-read.table(file = evall.buffer[i],header=T,row.names=1) - lista.buffer.final[[i]]<-lista.buffer[[i]]['TSS'] - lista.buffer.final[[i]][2]<-lista.buffer[[i]]['algoritmo'] + lista.buffer[[i]] <- read.table(file = evall.buffer[i], header = T, row.names = 1) + lista.buffer.final[[i]] <- lista.buffer[[i]]["TSS"] + lista.buffer.final[[i]][2] <- lista.buffer[[i]]["algoritmo"] } - mean.table.buffer<-matrix('NULL', nrow=length(lista.buffer.final)+2, ncol=dim(lista.buffer.final[[1]][1])[1]+2) - for (a in 1:length(lista.buffer.final)){ - for (b in 1:dim(lista.buffer.final[[1]][1])[1]){ - mean.table.buffer[a,b] <- as.matrix(lista.buffer.final[[a]][1])[b] + mean.table.buffer <- matrix("NULL", nrow = length(lista.buffer.final) + 2, ncol = dim(lista.buffer.final[[1]][1])[1] + + 2) + for (a in 1:length(lista.buffer.final)) { + for (b in 1:dim(lista.buffer.final[[1]][1])[1]) { + mean.table.buffer[a, b] <- as.matrix(lista.buffer.final[[a]][1])[b] } } - for (c in 1:length(mean.table.buffer[1,])){ - mean.table.buffer[nrow(mean.table.buffer)-1,c]<- mean(as.numeric((mean.table.buffer[,c])),na.rm=TRUE) - mean.table.buffer[nrow(mean.table.buffer),c]<- sd(as.numeric((mean.table.buffer[,c])),na.rm=TRUE) + for (c in 1:length(mean.table.buffer[1, ])) { + mean.table.buffer[nrow(mean.table.buffer) - 1, c] <- mean(as.numeric((mean.table.buffer[, + c])), na.rm = TRUE) + mean.table.buffer[nrow(mean.table.buffer), c] <- sd(as.numeric((mean.table.buffer[, + c])), na.rm = TRUE) } - colnames(mean.table.buffer)<- c(as.list(paste(lista.buffer.final[[1]][,2])),'stat','tratamento') - mean.table.buffer[,ncol(mean.table.buffer)-1]<-c(paste0('Com Buffer ',1:(nrow(mean.table.buffer)-2)),'mean','sd') - mean.table.buffer[,ncol(mean.table.buffer)]<-rep('Com Buffer') + colnames(mean.table.buffer) <- c(as.list(paste(lista.buffer.final[[1]][, 2])), + "stat", "tratamento") + mean.table.buffer[, ncol(mean.table.buffer) - 1] <- c(paste0("Com Buffer ", 1:(nrow(mean.table.buffer) - + 2)), "mean", "sd") + mean.table.buffer[, ncol(mean.table.buffer)] <- rep("Com Buffer") - #Resultados sem buffer---- - lista<-list() - lista.final<-list() + # Resultados sem buffer---- + lista <- list() + lista.final <- list() for (i in 1:length(evall)) { - lista[[i]]<-read.table(file = evall[i],header=T,row.names=1) - lista.final[[i]]<-lista[[i]]['TSS'] - lista.final[[i]][2]<-lista[[i]]['algoritmo'] + lista[[i]] <- read.table(file = evall[i], header = T, row.names = 1) + lista.final[[i]] <- lista[[i]]["TSS"] + lista.final[[i]][2] <- lista[[i]]["algoritmo"] } - mean.table<-matrix('NULL', nrow=length(lista.final)+2, ncol=dim(lista[[1]][1])[1]+2) - for (a in 1:length(lista.final)){ - for (b in 1:dim(lista.final[[1]][1])[1]){ - mean.table[a,b] <- as.matrix(lista.final[[a]][1])[b] + mean.table <- matrix("NULL", nrow = length(lista.final) + 2, ncol = dim(lista[[1]][1])[1] + + 2) + for (a in 1:length(lista.final)) { + for (b in 1:dim(lista.final[[1]][1])[1]) { + mean.table[a, b] <- as.matrix(lista.final[[a]][1])[b] } } - colnames(mean.table)<-c(as.list(paste(lista.final[[1]][,2])), 'stat','tratamento') - for (c in 1:length(mean.table[1,])){ - mean.table[nrow(mean.table)-1,c]<- mean(as.numeric((mean.table[,c])),na.rm=TRUE) - mean.table[nrow(mean.table),c]<- sd(as.numeric((mean.table[,c])),na.rm=TRUE) + colnames(mean.table) <- c(as.list(paste(lista.final[[1]][, 2])), "stat", "tratamento") + for (c in 1:length(mean.table[1, ])) { + mean.table[nrow(mean.table) - 1, c] <- mean(as.numeric((mean.table[, c])), + na.rm = TRUE) + mean.table[nrow(mean.table), c] <- sd(as.numeric((mean.table[, c])), na.rm = TRUE) } - mean.table[,ncol(mean.table)-1]<-c(paste0('sem Buffer ',1:(nrow(mean.table)-2)),'mean','sd') - mean.table[,ncol(mean.table)]<-rep('sem Buffer') + mean.table[, ncol(mean.table) - 1] <- c(paste0("sem Buffer ", 1:(nrow(mean.table) - + 2)), "mean", "sd") + mean.table[, ncol(mean.table)] <- rep("sem Buffer") - # Criando data.frame final ---- - #Finalizando mean - intermed_mean <- rbind(mean.table[nrow(mean.table)-1,], mean.table.buffer[nrow(mean.table.buffer)-1,]) + # Criando data.frame final ---- Finalizando mean + intermed_mean <- rbind(mean.table[nrow(mean.table) - 1, ], mean.table.buffer[nrow(mean.table.buffer) - + 1, ]) final_mean <- data.frame() - for (d in 1:(ncol(intermed_mean)-2)){ - #Sem bufer---- - final_mean[d,1]<-intermed_mean[1,d] - final_mean[d,2]<- colnames(intermed_mean)[d] - final_mean[d,3]<- intermed_mean[1,(ncol(intermed_mean))] - #com buffer ----- - final_mean[(ncol(intermed_mean)-2)+d,1]<-intermed_mean[2,d] - final_mean[(ncol(intermed_mean)-2)+d,2]<- colnames(intermed_mean)[d] - final_mean[(ncol(intermed_mean)-2)+d,3]<- intermed_mean[2,(ncol(intermed_mean))] + for (d in 1:(ncol(intermed_mean) - 2)) { + # Sem bufer---- + final_mean[d, 1] <- intermed_mean[1, d] + final_mean[d, 2] <- colnames(intermed_mean)[d] + final_mean[d, 3] <- intermed_mean[1, (ncol(intermed_mean))] + # com buffer ----- + final_mean[(ncol(intermed_mean) - 2) + d, 1] <- intermed_mean[2, d] + final_mean[(ncol(intermed_mean) - 2) + d, 2] <- colnames(intermed_mean)[d] + final_mean[(ncol(intermed_mean) - 2) + d, 3] <- intermed_mean[2, (ncol(intermed_mean))] } - #Finalizando SD - intermed_sd <- rbind(mean.table[nrow(mean.table),], mean.table.buffer[nrow(mean.table.buffer),]) + # Finalizando SD + intermed_sd <- rbind(mean.table[nrow(mean.table), ], mean.table.buffer[nrow(mean.table.buffer), + ]) final_sd <- data.frame() - for (d in 1:(ncol(intermed_sd)-2)){ - #Sem bufer---- - final_sd[d,1]<-intermed_sd[1,d] - final_sd[d,2]<- colnames(intermed_sd)[d] - final_sd[d,3]<- intermed_sd[1,(ncol(intermed_sd))] - #com buffer ----- - final_sd[(ncol(intermed_sd)-2)+d,1]<-intermed_sd[2,d] - final_sd[(ncol(intermed_sd)-2)+d,2]<- colnames(intermed_sd)[d] - final_sd[(ncol(intermed_sd)-2)+d,3]<- intermed_mean[2,(ncol(intermed_sd))] + for (d in 1:(ncol(intermed_sd) - 2)) { + # Sem bufer---- + final_sd[d, 1] <- intermed_sd[1, d] + final_sd[d, 2] <- colnames(intermed_sd)[d] + final_sd[d, 3] <- intermed_sd[1, (ncol(intermed_sd))] + # com buffer ----- + final_sd[(ncol(intermed_sd) - 2) + d, 1] <- intermed_sd[2, d] + final_sd[(ncol(intermed_sd) - 2) + d, 2] <- colnames(intermed_sd)[d] + final_sd[(ncol(intermed_sd) - 2) + d, 3] <- intermed_mean[2, (ncol(intermed_sd))] } - colnames(final_mean)<-c('mean','Algorithm','Tratamento') - colnames(final_sd)<-c('mean','Algorithm','Tratamento') + colnames(final_mean) <- c("mean", "Algorithm", "Tratamento") + colnames(final_sd) <- c("mean", "Algorithm", "Tratamento") - #removendo dados já utilizados: - #Buffer + # removendo dados já utilizados: Buffer rm(evall.buffer) rm(lista.buffer) rm(lista.buffer.final) rm(mean.table.buffer) - #s/buffer + # s/buffer rm(evall) rm(lista) rm(lista.final) rm(mean.table) - #demais + # demais rm(intermed_mean) rm(intermed_sd) -#Salvando em HD - if (file.exists(paste0(output.folder))==FALSE) dir.create(paste0(output.folder)) + # Salvando em HD + if (file.exists(paste0(output.folder)) == FALSE) + dir.create(paste0(output.folder)) - write.table(final_mean,file=paste0(output.folder,'/',sp,'_mean_results_tss.csv'), col.names=TRUE,row.names=FALSE) + write.table(final_mean, file = paste0(output.folder, "/", sp, "_mean_results_tss.csv"), + col.names = TRUE, row.names = FALSE) - write.table(final_sd,file=paste0(output.folder,'/',sp,'_sd_results_tss.csv'), col.names=TRUE,row.names=FALSE) - #Nao remover + write.table(final_sd, file = paste0(output.folder, "/", sp, "_sd_results_tss.csv"), + col.names = TRUE, row.names = FALSE) + # Nao remover rm(final_mean) rm(final_sd) rm(input.folder1) @@ -119,120 +131,101 @@ library("data.table") } #####----- -#Function to analyse(compare) the evaluate output from SDM, producing TSS Boxplot for each algorithm -consolidate.data <- function( - input.folder1="./TSS_Evaluate", - stat='mean', - name='TSS_analysis' -){ -library("data.table") - evall.list <- list.files(path = input.folder1, pattern=stat,full.names = T) #listing all evaluate files previusly organized. see read.eval() - lista<-list() +# Function to analyse(compare) the evaluate output from SDM, producing TSS +# Boxplot for each algorithm +consolidate.data <- function(input.folder1 = "./TSS_Evaluate", stat = "mean", name = "TSS_analysis") { + library("data.table") + evall.list <- list.files(path = input.folder1, pattern = stat, full.names = T) #listing all evaluate files previusly organized. see read.eval() + lista <- list() for (i in 1:length(evall.list)) { - lista[[i]]<-read.table(file = evall.list[i],header=T) + lista[[i]] <- read.table(file = evall.list[i], header = T) } - stats<-rbindlist(lista) #combining list - stats<-as.data.frame(stats) - #return(stats) - #Saving the boxplot graph in PNG + stats <- rbindlist(lista) #combining list + stats <- as.data.frame(stats) + # return(stats) Saving the boxplot graph in PNG library(ggplot2) - ggplot(stats, aes(x=Tratamento, y=mean, fill=Tratamento)) + geom_boxplot(notch=TRUE) + facet_grid(~Algorithm) + - labs(y='') + - theme(axis.title.x=element_blank(), - axis.text.x=element_blank(), - axis.ticks.x=element_blank(), - axis.text.y=element_text(size=rel(1.5)), - strip.text = element_text(face="bold", size=rel(1.25)), - legend.title=element_text(face='bold',size=rel(1.5)), - legend.text=element_text(size=rel(1.5)), - legend.position="bottom") - ggsave(filename=paste0(input.folder1,"/",name,".png")) -# dev.off() #Not necessary when using ggsave -#return(stats) + ggplot(stats, aes(x = Tratamento, y = mean, fill = Tratamento)) + geom_boxplot(notch = TRUE) + + facet_grid(~Algorithm) + labs(y = "") + theme(axis.title.x = element_blank(), + axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.text.y = element_text(size = rel(1.5)), + strip.text = element_text(face = "bold", size = rel(1.25)), legend.title = element_text(face = "bold", + size = rel(1.5)), legend.text = element_text(size = rel(1.5)), legend.position = "bottom") + ggsave(filename = paste0(input.folder1, "/", name, ".png")) + # dev.off() #Not necessary when using ggsave return(stats) rm(stats) rm(lista) rm(evall.list) } #####----- -#Function to save the species name and final occurrence. -species.table <- function( - output.folder="./TSS_Evaluate", - table.name="n.sp.modelado"){ -n.sp<-table(occs$sp) -n.sp<-as.data.frame(n.sp) -#head(n.sp) -#dim(n.sp) -n.sp <- subset(n.sp, Freq > 10) -#dim(n.sp) -n.sp.modelado <- subset(n.sp, Var1 %in% (names.sp[N])) -#head(n.sp.modelado) -#dim(n.sp.modelado) -write.csv(n.sp.modelado, paste0(output.folder,'/',table.name,'.csv')) -rm(n.sp) -rm(n.sp.modelado) +# Function to save the species name and final occurrence. +species.table <- function(output.folder = "./TSS_Evaluate", table.name = "n.sp.modelado") { + n.sp <- table(occs$sp) + n.sp <- as.data.frame(n.sp) + # head(n.sp) dim(n.sp) + n.sp <- subset(n.sp, Freq > 10) + # dim(n.sp) + n.sp.modelado <- subset(n.sp, Var1 %in% (names.sp[N])) + # head(n.sp.modelado) dim(n.sp.modelado) + write.csv(n.sp.modelado, paste0(output.folder, "/", table.name, ".csv")) + rm(n.sp) + rm(n.sp.modelado) } #####----- -#Functio for making and cropping rasterLayer -rasterCrop<- function(sp, - input.folder = './buffermax', - mascara = raster, - crop = vctor_shape -){ +# Functio for making and cropping rasterLayer +rasterCrop <- function(sp, input.folder = "./buffermax", mascara = raster, crop = vctor_shape) { library(raster) - rasters <- - list.files( - path = paste0(input.folder,"/",sp), pattern = ".tif$",full.names = T - ) + rasters <- list.files(path = paste0(input.folder, "/", sp), pattern = ".tif$", + full.names = T) mod <- NULL for (i in 1:length(rasters)) { - ?raster + `?`(raster) mod <- raster(rasters[i]) - #plot(mod) + # plot(mod) mod <- mask(mod, mascara) mod <- crop(mod, crop) - writeRaster(mod, filename=rasters[i], overwrite = T,format = "GTiff") + writeRaster(mod, filename = rasters[i], overwrite = T, format = "GTiff") rm(mod) - }} + } +} -#Function to analyse perc ensemble between partitions and algorithms---- -partitionEnsemble <- function( - sp, - projeto="./FLORA_buffermax2", - input.folder='presfinal', - algoritmos = c("maxent", "rf", "svm")){ +# Function to analyse perc ensemble between partitions and algorithms---- +partitionEnsemble <- function(sp, projeto = "./FLORA_buffermax2", input.folder = "presfinal", + algoritmos = c("maxent", "rf", "svm")) { library(raster) library(data.table) - for (algo in 1:length(algoritmos)){ + for (algo in 1:length(algoritmos)) { - #Lista raster binario saida 7 do algoritmo - modelos.bin <- list.files(path = paste0(projeto,"/",sp,"/", input.folder), full.names = T,pattern = paste0("bin7",sp,algoritmos[algo],'.tif')) - if(length(modelos.bin)==0){ + # Lista raster binario saida 7 do algoritmo + modelos.bin <- list.files(path = paste0(projeto, "/", sp, "/", input.folder), + full.names = T, pattern = paste0("bin7", sp, algoritmos[algo], ".tif")) + if (length(modelos.bin) == 0) { - ensemble<-data.frame(ensemble=NA) - ensemble$n.sp<-table2[(which(table2$sp==sp)),2] + ensemble <- data.frame(ensemble = NA) + ensemble$n.sp <- table2[(which(table2$sp == sp)), 2] ensemble$algoritmo <- algoritmos[algo] ensemble$nome.sp <- sp ensemble_list[[algo]] <- ensemble - } else{ + } else { - #carrega raster + # carrega raster mod.bin <- raster(modelos.bin) - #estatistica quantidade de pixels - zonal.result <- zonal(zonal_maks, mod.bin, fun='count', na.rm=TRUE, digits=2) + # estatistica quantidade de pixels + zonal.result <- zonal(zonal_maks, mod.bin, fun = "count", na.rm = TRUE, + digits = 2) zonal.result <- as.data.frame(zonal.result) - #separa apenas resultado impoirtante - ensemble <- zonal.result[nrow(zonal.result),ncol(zonal.result)]/sum(zonal.result[c(2:nrow(zonal.result)),ncol(zonal.result)]) + # separa apenas resultado impoirtante + ensemble <- zonal.result[nrow(zonal.result), ncol(zonal.result)]/sum(zonal.result[c(2:nrow(zonal.result)), + ncol(zonal.result)]) rm(zonal.result) - #Construindo resultados - ensemble<-as.data.frame(ensemble) - ensemble$n.sp<-table2[(which(table2$sp==sp)),2] + # Construindo resultados + ensemble <- as.data.frame(ensemble) + ensemble$n.sp <- table2[(which(table2$sp == sp)), 2] ensemble$algoritmo <- algoritmos[algo] ensemble$nome.sp <- sp ensemble_list[[algo]] <- ensemble @@ -243,48 +236,48 @@ partitionEnsemble <- function( return(ensemble_list) } -#Function to analyse the percentual ensemble between <> algorithms/species---- -spEnsemble <- function( - sp, - projeto="./FLORA_buffermax2", - input.folder='ensemble'){ +# Function to analyse the percentual ensemble between <> algorithms/species---- +spEnsemble <- function(sp, projeto = "./FLORA_buffermax2", input.folder = "ensemble") { library(raster) library(data.table) - #Lista raster binario saida 7 do algoritmo - modelos.bin <- list.files(path = paste0(projeto,"/",sp,"/", input.folder), full.names = T,pattern = paste0(".bin7",'_ensemble.tif')) - if(length(modelos.bin)==0){ - ensemble<-data.frame(ensemble=NA) - ensemble$n.sp<-table2[(which(table2$sp==sp)),2] + # Lista raster binario saida 7 do algoritmo + modelos.bin <- list.files(path = paste0(projeto, "/", sp, "/", input.folder), + full.names = T, pattern = paste0(".bin7", "_ensemble.tif")) + if (length(modelos.bin) == 0) { + ensemble <- data.frame(ensemble = NA) + ensemble$n.sp <- table2[(which(table2$sp == sp)), 2] ensemble$sp <- sp - }else{ + } else { - #carrega raster + # carrega raster mod.bin <- raster(modelos.bin) - #estatistica quantidade de pixels - zonal.result <- zonal(zonal_maks, mod.bin, fun='count', na.rm=TRUE, digits=2) + # estatistica quantidade de pixels + zonal.result <- zonal(zonal_maks, mod.bin, fun = "count", na.rm = TRUE, digits = 2) zonal.result <- as.data.frame(zonal.result) - #separa apenas resultado impoirtante - ensemble <- zonal.result[nrow(zonal.result),ncol(zonal.result)]/sum(zonal.result[c(2:nrow(zonal.result)),ncol(zonal.result)]) + # separa apenas resultado impoirtante + ensemble <- zonal.result[nrow(zonal.result), ncol(zonal.result)]/sum(zonal.result[c(2:nrow(zonal.result)), + ncol(zonal.result)]) rm(zonal.result) - #Construindo resultados - ensemble<-as.data.frame(ensemble) - ensemble$n.sp<-table2[(which(table2$sp==sp)),2] + # Construindo resultados + ensemble <- as.data.frame(ensemble) + ensemble$n.sp <- table2[(which(table2$sp == sp)), 2] ensemble$sp <- sp } return(na.omit(ensemble)) } # Funtion to conver raaster to smaller type---- -rast.convert<-function(sp, - projeto='./FLORA_buffermax2', - input.folder="ensemble"){ - modelos.bin <- list.files(path = paste0(projeto,"/",sp,"/", input.folder), full.names = T,pattern = paste0(".bin7_ensemble50.tif")) - if(length(modelos.bin)==0){} else { - modelos.bin<-raster(modelos.bin) - writeRaster(modelos.bin, filename = paste0("./",projeto,"/final/",sp,".tif"), datatype="INT1U") +rast.convert <- function(sp, projeto = "./FLORA_buffermax2", input.folder = "ensemble") { + modelos.bin <- list.files(path = paste0(projeto, "/", sp, "/", input.folder), + full.names = T, pattern = paste0(".bin7_ensemble50.tif")) + if (length(modelos.bin) == 0) { + } else { + modelos.bin <- raster(modelos.bin) + writeRaster(modelos.bin, filename = paste0("./", projeto, "/final/", sp, + ".tif"), datatype = "INT1U") } -} \ No newline at end of file +} diff --git a/ENM/fct/write_Raster.R b/ENM/fct/write_Raster.R index c5e15f1..2c0d75d 100644 --- a/ENM/fct/write_Raster.R +++ b/ENM/fct/write_Raster.R @@ -1,7 +1,8 @@ -write_Raster <- function(x,mask=mask){ - if (class(mask) == "RasterLayer"){ - x <- mask(x , mask) - x <- crop(x , mask) - writeRaster(x,filename=paste0(output.folder,"/",sp,"/",proj,"/",algo,"_cont_",sp,"_",i,".tif"),overwrite=T,datatype="INT1U") +write_Raster <- function(x, mask = mask) { + if (class(mask) == "RasterLayer") { + x <- mask(x, mask) + x <- crop(x, mask) + writeRaster(x, filename = paste0(output.folder, "/", sp, "/", proj, "/", + algo, "_cont_", sp, "_", i, ".tif"), overwrite = T, datatype = "INT1U") } -} \ No newline at end of file +}