Modelagem de Distribuição de Espécies

Os dados básicos para a modelagem potencial de espécies são:

• Dados bióticos
• Dados abióticos

Para facilitar, utilizaremos dados bióticos (bradypus.csv) e abióticos existentes no pacote DISMO.

Dados Bióticos:

file <- paste(system.file(package="dismo"), "/ex/bradypus.csv", sep="")
bradypus <- read.table(file, header=TRUE, sep= ",")
head(bradypus)

              species    lon    lat
1 Bradypus variegatus -65.40 -10.38
2 Bradypus variegatus -65.38 -10.38
3 Bradypus variegatus -65.13 -16.80
4 Bradypus variegatus -63.67 -17.45
5 Bradypus variegatus -63.85 -17.40
6 Bradypus variegatus -64.42 -16.00

#Deixaremos o objeto 'bradypus' apenas com as coordenadas
bradypus <- bradypus[,2:3]

Dados abióticos

##dados abióticos (ambientais)
#criamos o caminho para os rasters
files2 <- list.files(path=paste(system.file(package="dismo"),
'/ex', sep=''), pattern='grd', full.names=TRUE )
#Definimos os rasters como objeto 'predictors'
predictors <- stack(files2)

Uma vez que temos como camada ambiental a informação de bioma e esta não será utilizada na modelagem, vamos desconsidera-la.

pred_nf <- dropLayer(predictors,'biome')

Agora que temos os dados bióticos e abióticos, precisamos extrair os valores ambientais das localidades onde há registros de ocorrência, usando a função 'extract'.
A tais valores daremos o nome de 'presvals':

presvals <- extract(pred_nf, bradypus)

Faremos então a divisão dos nossos dados de bradypus em 5 grupos diferentes para utilizarmos 4 para treino e apenas um para teste

group <- kfold(bradypus, 5)
pres_train <- bradypus[group != 1, ]
pres_test <- bradypus[group == 1, ]

Agora criaremos pontos leatórios, como bakcground e os dividiremos em 5 grupos diferentes: 4 para treino e um para teste:

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

Rodando a modelagem:

Com algoritmo Bioclim

bc <- bioclim(pred_nf, pres_train)

Com algoritmo Mahalanobis

mm <- mahal(pred_nf, 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, pred_nf) 
e

class          : ModelEvaluation 
n presences    : 23 
n absences     : 200 
AUC            : 0.7717 
cor            : 0.2312 
max TPR+TNR at : 0.08592 

Avaliando o modelo gerado pelo Mahalanobis:

e_mm <- evaluate(pres_test, backg_test, mm, pred_nf)
e_mm

class          : ModelEvaluation 
n presences    : 23 
n absences     : 200 
AUC            : 0.8633 
cor            : 0.1223 
max TPR+TNR at : -3.349 

Definindo um threshold: Boiclim

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

[1] 0.08592

Voltando para o espaço geográfico: Bioclim:

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

class       : RasterLayer 
dimensions  : 112, 116, 12992  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -90, -32, -33, 23  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.7742  (min, max)

Voltando para o espaço geográfico:
Mahalanobis:

pm = predict(pred_nf, mm, ext=ext, progress='')
pm[pm < -10] <- -10
pm

class       : RasterLayer 
dimensions  : 112, 116, 12992  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -90, -32, -33, 23  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : -10, 1  (min, max)

Definindo um threshold: Mahalanobis

threshold_mm <- threshold(e_mm,'spec_sens')
threshold_mm

[1] -3.349

Modelos binários definidos pelos threshold: