Curso Atualidades em Macroecologia
Escola de Altos Estudos CAPES - Brasil
Universidade Federal de Sergipe
dismo (Hijmans & Elith 2011) e código do meu amigo Vijay Barve (https://vijaybarve.wordpress.com) para mapear com ggplot2 e Google MapsA gente já viu de onde vem os dados primários de biodiversidade (e.g. dados de ocorrência; armazenados em museus, herbaria, coleções científicas). Muitos desses dados estão agora disponíveis na internet!!!
Carregar os pacotes (TODOS!) necessarios pra esse exercicio (eles já devem estar instalados nos seus computadores)
library(sp)
library(raster)
library(maptools)
library(maps)
library(dismo)
library(rgdal)
library(ggmap)
library(ggplot2)
library(rgbif)Obter dados de ocorrência de uma espécie qualquer a partir de base da dados GBIF (data.gbif.org) direitamente desde o software R Nota: tem que estar conectados à internet. Pode demorar dependendo se a espécie tem (ou não) muitos registros
sp_1 <- gbif("Chrysocyon","brachyurus")Nota: vcs podem colocar o nome do genero e a espécie desejada.
Olhar as dimensões desse objeto (i.e. quantas linhas e quantas colunas)
dim(sp_1)Conferir os nomes das colunas (pra depois procurar só a longitude e a latitude)
names(sp_1)Crear outro objeto a partir do anterior pra ficar com long/lat apenas
sp1_points <- sp_1[,c("lon","lat")]Nota: o nome da variavel pode ser diferente (.e.g “LATITUDE”, “Latidude”, “lat”, etc. Sempre olhar antes)
Graficar esses pontos (x,y)
plot(sp1_points[,1],sp1_points[,2],col="blue",pch=19)Adicionar o mapa do mundo. Pra isso, precisam carregar esse mapa (está no pacote “maptools”)
data(wrld_simpl)
dev.new()
plot(wrld_simpl,add=T)Tem algum dado estranho??
Graficar primeiro o mundo e depois os pontos
plot(wrld_simpl,col="light yellow")Adicionar os pontos
points(sp1_points[,1],sp1_points[,2],col="blue",cex=0.2)Nota: os argumetos “pch” acima e “cex” aqui, apenas definem o tamanho e tipo dos pontinhos. Podem brincar com isso e achar o melhor
Dependendo da espécies, o mapa pode ser muito grande né? Da pra fazer um menor sabendo onde que os pontos estão. Nesse caso, vamos crear um mapa só pra o Brasil (ou pra qualquer pais, só mudar o “NAME”)
brazil_map <- wrld_simpl[wrld_simpl$NAME=="Brazil",]Agora podem repetir os passos anteriores e mapear apenas pra o Brasil.
conferir se tem “duplicatas”
sp1_dups <- duplicated(sp1_points)
### NOTA: a função "duplicated" retorna o resultado de um teste logico (e.g. TRUE or FALSE)
# quantos são duplicatas?
length(which(sp1_dups==TRUE))
# quantos não são duplicatas?
length(which(sp1_dups==FALSE))
# manter o as linhas que não são duplicatas
sp1_dups_row <- which(sp1_dups==TRUE)
# gerar outro "objeto" sem as duplicatas
sp1_nodups <- sp1_points[-sp1_dups_row,]
# quais são as dimensões do novo objeto?
dim(sp1_nodups)
# olhar os primeiros elementos/dados do novo objeto
head(sp1_nodups)conferir se tem dados com longitude/latitude zero
sp1_lonzero = subset(sp1_nodups, lon==0)
#conferir se tem dados com latitude zero
sp1_latzero = subset(sp1_nodups, lat==0)Nota: lembrar que os nomes podem ser diferentes (Latitude, latitude, lat, etc…)
conferir se tem dados com NA’s
sp1_nas <- which(is.na(sp1_nodups[,1])==TRUE)
sp1_points_nonas <- sp1_nodups[-sp1_nas,]E agora? Ficou igual/melhor? A gente sabe que esse bicho não está em alguns dos lugares aonde aparecem pontos, certo? Vamos tira-los
sp1_points_final <- sp1_points_nonas[-which(sp1_points_nonas[,2]>0),]Baixar dados pra uma outra espécie
sp_2 <- gbif("Leopardus","colocolo")Ficar apenas com long/lat (lembrem que tem que olhar os “names” do objeto sp_2 pra saber quais colunas tem essa info)
sp2_points <- sp_2[,c("lon","lat")]Graficar
plot(wrld_simpl, col = "light yellow", axes = F)
points(sp2_points[,1], sp2_points[,2], col = "red", cex = 0.5)
text(-140, -50, "Leopardus\ colocolo")Graficar com outros mapas, usando outros dados (e.g. baixar pra outra espécie usando outro pacote “rgbif”).
jaguarundi <- occ_data(scientificName = 'Herpailurus yagouaroundi', hasCoordinate = TRUE)Nota: vcs podem pegar o dado/especie que quiser.
Obter o mapa do mundo
world = map_data("world")Graficar os pontos usando o mapa como poligono e definendo caracteristicas a esse ponto (color, tamanho, etc)
ggplot(world, aes(long, lat)) +
geom_polygon(aes(group = group), fill = "white",
color = "gray40", size = .2) +
geom_jitter(data = jaguarundi$data,
aes(jaguarundi$data$decimalLongitude, jaguarundi$data$decimalLatitude), alpha=0.6,
size = 4, color = "red") +
labs(title = "Herpailurus yagouarundi")Ficou legal? e usando outra espécie?
Agora com mapas do GoogleMaps
aardvark = occ_data(scientificName = 'Orycteropus afer', hasCoordinate = TRUE)
# definir o "extent", a "janela" do mapa
e = extent( -179 , 179 , -80 , 80 )
# pegar essa janela do mapa desde o GoogleMaps
r = gmap(e)
# graficar
plot(r, interpolate=TRUE, main="Map")
# juntar os valores de long/lat
xy1=cbind(aardvark$data$decimalLongitude,aardvark$data$decimalLatitude)
# adicionar os pontos, usando a mesma projeção que o mapa (i.e. "Mercator")
points(Mercator(xy1) , col='blue', pch=20)
# colocar o nome da especie
text(160,0, "\n\n\nOrycteropus \nafer", adj = c(0,1),
col="red")Fazer “zoom” com esse dados e mapear de novo
wmap1 = qmap('Africa',zoom=2)
wmap1 +
geom_jitter(data = aardvark$data,
aes(decimalLongitude, decimalLatitude),
alpha=0.6, size = 4, color = "red") +
labs(title = "Orycteropus afer")Essa parte aquí pode ser usada pra gerar mapas “simples” baseados na geometría (e.g. distancias entre pontos, etc.). Não levam em consideração as variáveis ambientais.
This model draws a convex hull around all ’presence’ points.
#crear um raster que vai ser usado como "máscara"
r <- raster()
values(r) <- 1
ext <- c(-100,-30,-60,14)
r1 <- crop(r,ext)
guara_hull <- convHull(sp1_points_final)
guara_hullmodel <- predict(guara_hull, r1)
plot(guara_hullmodel, main='Convex Hull')
plot(wrld_simpl, add=TRUE, border='dark grey')
points(sp1_points_final, pch='+')This model draws circles around all ’presence’ points.
guara_circ <- circles(sp1_points_final, lonlat=TRUE)
guara_circles <- predict(r1, guara_circ, mask=TRUE)
plot(guara_circles, main='Circles')
points(sp1_points_final, pch=19, cex=0.5, col="red")
#"brincar" com a distância usada pra desenhar os "circles"
guara_circ2 <- circles(sp1_points_final, d=1000000, lonlat = TRUE)
guara_circles2 <- predict(r1, guara_circ2, mask=TRUE)
plot(guara_circles2, main='Circles')
points(sp1_points_final, pch=19, cex=0.5, col="red")Agora, vamos pegar dados ambientais/climáticos
Baixar os dados de clima do site www.worldclim.org (pode demorar um pouquinho, par de minutos ou menos)
world_bioclim<-getData("worldclim", var = "bio", res = 10)Essas variaveis estão todas “juntas”. Vamos separar e deixar em um objeto aparte que pode ser usado pra pegar uma por uma
world_bioclim_inds<-unstack(world_bioclim)BIO1 = Annual Mean Temperature
BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
BIO3 = Isothermality (BIO2/BIO7) (* 100)
BIO4 = Temperature Seasonality (standard deviation *100)
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO7 = Temperature Annual Range (BIO5-BIO6)
BIO8 = Mean Temperature of Wettest Quarter
BIO9 = Mean Temperature of Driest Quarter
BIO10 = Mean Temperature of Warmest Quarter
BIO11 = Mean Temperature of Coldest Quarter
BIO12 = Annual Precipitation
BIO13 = Precipitation of Wettest Month
BIO14 = Precipitation of Driest Month
BIO15 = Precipitation Seasonality (Coefficient of Variation)
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
BIO18 = Precipitation of Warmest Quarter
BIO19 = Precipitation of Coldest Quartergraficar uma das variaveis
plot(world_bioclim_inds[[1]], zlim=c(-270,320), axes = F, box=F, col=rev(heat.colors(25)), ext = matrix(c(180,90,-180,-90), nrow = 2))O um plot mais “chick”
library(rasterVis)
gplot(world_bioclim_inds[[1]],maxpixels=100000)+ # reference the data
geom_raster(aes(fill=value)) + scale_fill_gradientn(colours=rev(c("brown","red","yellow","darkgreen","green")))Tem que graficar alguma dessas variaveis apenas pra America e pra o Brasil. Pra isso tem que mexer no “ext=matrix(c(XX,XX,XX,XX)))” da função acima, colocando as coordenadas certas pra America e o Brasil (Por separado, um mapa pra cada)
vamos ficar com duas variáveis apenas
bio1<- world_bioclim_inds[[1]]
bio12<- world_bioclim_inds[[12]]
env.variables<- stack (bio1, bio12)Estabelecer o “extent” da América do Sul e cortar as variáveris ambientais nesse extent
am.extent <- c(-100,-30,-60,14)
am.env <- crop(env.variables,am.extent)Gerar o modelo Bioclim (Ecological Niche Model: ENM)
# Converter os pontos da espécie em uma `matrix`
sp1_points_mat <- as.matrix(sp1_points_final)
# Rodar o modelo ENM (Bioclim)
sp1_bioclim <- bioclim(am.env, sp1_points_mat)
plot(sp1_bioclim, a=1, b=2, p=0.85)
# Gerar uma predição da área adequada ("suitable") baseada no ENM
sp1.map<- predict (am.env, sp1_bioclim)
# Plotar o resultado da predição baseada no ENM
plot (sp1.map)
# Avaliar o modelo ENM
group <- kfold(sp1_points_final, 5)
pres_train <- sp1_points_final[group != 1, ]
pres_test <- sp1_points_final[group == 1, ]
backg <- randomPoints(am.env, n=1000, ext=am.extent, extf = 1.25)
colnames(backg) = c('lon', 'lat')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]
e <- evaluate(pres_test, backg_test, sp1_bioclim, am.env)
# Estabelecer um "threshold" pra cortar a predição e gerar um mapa binario
# max TPR+TNR: maximum of the sum of the sensitivity (true positive rate) and specificity (true negative rate). This is sometimes used as a threshold for setting cells to presence or absence.
threshold <- e@t[which.max(e@TPR + e@TNR)]
plot(sp1.map > threshold, main='presence/absence')
# Adicionar os pontos reais da espécie
points(sp1_points_mat[,1],sp1_points_mat[,2],cex=0.2)Como se ve?
Tem regiões aonde a espécie foi predita mas ta ausente?
NOTA: Salvar o “workspace” do R. Isto é tudo o que foi feito neste exercicio tem que dicar guardado como um arquivo “.Rdata”. Esses objetos vão ser usado em outros exercicios. NÃO FECHEM O R SEM SALVAR O ARQUIVO!!!