Modelagem de Distribuição de Espécies

Trabalhando com dados próprios

Professora: Marinez Ferreira de Siqueira
Andrea Sánchez-Tapia
Pablo Viany Prieto
Felipe Sodré Barros
Rafael Oliveira
2014

Primeiro passo: Habilitar os pacotes a serem usados:

library(dismo)
library(raster)

Para usarmos os dados bióticos de cada alunos, identificaremos o workingDirectiory do R onde a tabelas salva em formato .csv ficará armazenada:

getwd()

[1] "/home/iis/Dropbox/SDM2"

##Caso seja necessário, defina o wd
#setwd(u:/aluno/)

Acessando e lendo dados bióticos:

dados_bio <- read.table("caryocar.csv", header=TRUE, sep= ",")
head(dados_bio)

               especie   long    lat
1 Caryocar brasiliense -49.25 -20.55
2 Caryocar brasiliense -47.52 -14.12
3 Caryocar brasiliense -46.25 -20.92
4 Caryocar brasiliense -48.92 -10.08
5 Caryocar brasiliense -51.77 -12.82
6 Caryocar brasiliense -48.96 -22.32

#Deixaremos o objeto importado em csv apenas com as coordenadas.
#Considerando que Lat/Longocupam a 2 e 3 coluna...
dados_modelo <- dados_bio[,2:3]

Dados abióticos

#criamos o caminho para os rasters
rasters <- list.files(path=paste(getwd(),'/datasets', sep=''), pattern='bio', full.names=TRUE )
#Definimos os rasters como objeto 'predictors'
predictors <- stack(rasters)

plot(predictors)

plot of chunk unnamed-chunk-5

Vamos visualizar os dados bióticos espacializados sobre uma das camadas preditoras:

plot(predictors, 1)
points(dados_modelo, bg='red', cex=1,pch=21)

plot of chunk unnamed-chunk-6

Extraindo os valores ambientais das localidades onde há registros de ocorrência, usando a função 'extract':

presvals <- extract(predictors, dados_modelo)

Como usaremos o algoritmo Bioclim, que não precisa de dados de ausência, não geraremos random points nesta estapa.

#ext <- extent(-90, -32, -33, 23)
#backgr <- randomPoints(predictors, 500)
#absvals <- extract(pred_nf, backgr)

Antes de realizar a modelagem, precisamos definir o como tais dados serão particionadas para treino e teste, de presença apenas.

Ou seja, particionaremos os dados de presença em dois: Uma parte (treino) será usada para treinar o algoritmos e realizar a modelagem. A outra partição (teste) será usada para testar se a modelagem possui erros de comissão e/ou erros de omissão.

group <- kfold(dados_modelo, 3)
pres_train <- dados_modelo[group != 1, ]
pres_test <- dados_modelo[group == 1, ]

Agora a divisão dos nossos dados em 3 grupos diferentes para utilizarmos 2 para treino e apenas um para teste:

ext <- extent(-184, 184, -60, 90)
backg <- randomPoints(predictors, n=1000, ext=ext, extf = 1.25)
colnames(backg) = c( 'lon' ,  'lat' )
##Criando e separando as partições:
group <- kfold(backg, 3)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]

Para visualizar onde estão as ocorrências dos seus dados e onde estão os pontos aleatórios:

r<- raster(predictors, 1)
plot(!is.na(r), col=c( 'white' ,  'light grey' ), legend=FALSE)
points(pres_train, pch=  '+' , col= 'green' )
points(pres_test, pch= '+' , col= 'blue' )

plot of chunk unnamed-chunk-11

Dados aleatórios

plot(!is.na(r), col=c( 'white' ,  'light grey' ), legend=FALSE)
points(backg_train, pch= '-' , cex=0.5, col= 'yellow' )
points(backg_test, pch= '-' ,  cex=0.5, col= 'black' )

plot of chunk unnamed-chunk-12

Rodando a modelagem com algoritmo Bioclim:

bc <- bioclim(predictors, pres_train)

Avaliando o modelo gerado pelo Bioclim, usando o conjunto de presença deixada para teste e os aleatórios de teste:

e <- evaluate(pres_test, backg_test, bc, predictors) 
e

class          : ModelEvaluation 
n presences    : 38 
n absences     : 333 
AUC            : 0.9397 
cor            : 0.6657 
max TPR+TNR at : 0.01306

Definindo um threshold: Boiclim

threshold <- e@t[which.max(e@TPR + e@TNR)]
threshold

[1] 0.01306

Voltando para o espaço geográfico: Bioclim

pb <- predict(predictors, bc, ext=ext, progress='')
pb

class       : RasterLayer 
dimensions  : 1470, 3675, 5402250  (nrow, ncol, ncell)
resolution  : 0.1, 0.1  (x, y)
extent      : -183.5, 184, -58.73, 88.27  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
data source : in memory
names       : layer 
values      : 0, 0.8684  (min, max)

Visualizando o resultado no espaço geográfico: Bioclim Modelo Originai:

plot(pb, main= 'Bioclim, raw values' )
points(pres_train, pch= '+' )

plot of chunk unnamed-chunk-17

Modelo binário:

plot(pb > threshold, main= 'Bioclim presence/absence ')
points(pres_train, pch= '+' )

plot of chunk unnamed-chunk-18

Exportando o modelo gerado:

pb_final<-(pb > threshold)
writeRaster(pb_final,filename='Bioclim_Threshold.asc', format='ascii', datatype='INT2U', overwrite=T)

class       : RasterLayer 
dimensions  : 1470, 3675, 5402250  (nrow, ncol, ncell)
resolution  : 0.1, 0.1  (x, y)
extent      : -183.5, 184, -58.73, 88.27  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/iis/Dropbox/SDM2/Bioclim_Threshold.asc 
names       : Bioclim_Threshold 
values      : 0, 1  (min, max)