Modelos de distribuição de espécies (SDM, do inglês Species Distribution Modelling) podem ser gerados no R com o pacote dismo (Hijmans et al. 2017). Para uma revisão extensa sobre SDM usando dismo, veja aqui.
Para fazer as predições sobre distribuição das espécies precisamos de camadas ambientais e dados de ocorrência das espécies desejadas. As camadas aqui utilizadas são do WorldClim (Fick & Hijmans 2017), de onde podem ser baixadas. Porém, podemos baixá-las usando a função getData() do pacote raster.
As ocorrências de espécies podem ser baixadas pelo GBIF ou também pelo R, usando a função gbif() do pacote dismo. As ocorrências podem conter duplicatas, problemas com coordenadas geográficas, etc. Por isso, vamos usar o pacote CoordinateCleaner (Zizka et al. 2019) que compara as coordenadas das ocorrências com várias fontes de informação, e.g. oceanos, centro de pesquisas, centroide dos países. Ou seja, se a ocorrência de uma planta estiver fora do continente (no oceano), no centróide de um país, num centro de pesquisa (ou jardim botânico por exemplo onde é cultivada), elas são apontadas e o usuário pode retirá-las dos dados.
Vamos começar carregando todos pacotes que utilizaremos nesta rotina:
> library(raster);library(rgdal);library(dismo);library(maps)
> library(maptools); library(CoordinateCleaner)
Em seguida, vamos baixar dados de ocorrência do Buriti para a América do Sul. Para tanto, criamos um raster com os limites do continente usando a função extent:
> ext<-extent(-81.5, -34, -57, 13) #America do Sul
Agora baixamos os dados criando o objeto Mauritia, em seguida, vamos eliminar duplicatas com a função subset, e por fim escolhemos só os registros que são "PRESERVED_SPECIMEN" ou seja, exclui observações não coletadas ou registros fósseis.
> Mauritia<-gbif("Mauritia flexuosa",download=T,geo=T, ext=ext)
> Mauritia<-subset(Mauritia, !duplicated(subset(Mauritia, select=c(lon,lat))))
> Mauritia<-subset(Mauritia, basisOfRecord=="PRESERVED_SPECIMEN")
Agora criamos uma nova tabela Mauritia com nomes de colunas diferentes (exigidas pelo CoordinateCleaner):
> Mauritia<-data.frame(species= Mauritia$acceptedScientificName, decimallongitude=Mauritia$lon,decimallatitude=Mauritia$lat)
> Mauritia<-subset(Mauritia, !decimallongitude=="NA") #caso alguma Lat ou Long exista, a linha é excluida
> flags<-clean_coordinates(x=Mauritia, tests=c("capitals","centroids", "equal", "gbif", "institutions", "outliers", "seas", "zeros"))
> Mauritia<-subset(flags,.summary==TRUE)
> Mauritia<-data.frame(Long= Mauritia$decimallongitude,Lat= Mauritia$decimallatitude)
Agora fazemos a mesma rotina acima para o gênero Mauritiella:
> Mauritiella<-gbif("Mauritiella",download=T,geo=T, ext=ext) #import data
> Mauritiella<-subset(Mauritiella, !duplicated(subset(Mauritiella, select=c(lon,lat))))
> Mauritiella<-subset(Mauritiella, basisOfRecord== "PRESERVED_SPECIMEN")
> Mauritiella<-data.frame(species= Mauritiella$acceptedScientificName, decimallongitude=Mauritiella$lon, decimallatitude=Mauritiella$lat)
> Mauritiella<-subset(Mauritiella, !decimallongitude=="NA")
> flags<-clean_coordinates(x=Mauritiella, tests=c("capitals","centroids", "equal", "gbif", "institutions", "outliers", "seas", "zeros")) #record-level tests
> Mauritiella<-subset(flags,.summary==TRUE)
> Mauritiella<-data.frame(Long= Mauritiella$decimallongitude,Lat= Mauritiella$decimallatitude)
Baixando os dados do WorldClim
Baixamos os dados do worldclim com a função getData. Ela baixa 19 camadas bioclimáticas que podem ser conferias aqui. A camada topográfica também podem ser baixada com o getData, mas é feita separadamente.
> r <- getData("worldclim",var="bio",res=10)
> alt <- getData('worldclim', var='alt', res=10)
Este procedimento cria uma pasta no seu diretório chamada wc10 (WorldClim, resolução 10) com os arquivos de cada camada. Os arquivos zipados deve ser excluidos. Para ler os dados para o console, podemos listar os arquivos da seguinte forma:
> rasters<-list.files(pattern=".bil")
Todas camadas são listadas no objeto rasters. Porém, vamos escolher apenas algumas camadas, fazemos isso indicando quais queremos:
> rasters<-c("alt.bil", "bio1.bil", "bio5.bil", "bio6.bil", "bio12.bil","bio16.bil","bio17.bil")
As camadas escolhidas são:
BIO1 = Annual Mean Temperature
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO12 = Annual Precipitation
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
Agora é necessário sobrepor todas camadas numa só, fazemos isso com a função raster:
> stackrasters<-raster(rasters[1])
> for(x in rasters){stackrasters<-stack(stackrasters,raster(x))}
Um problema ocorre neste passo, a primeira camada é repetida. Visualize a sobreposição de camadas, veja que há duas camadas "alt."
> stackrasters
class : RasterStack
dimensions : 900, 2160, 1944000, 8 (nrow, ncol, ncell, nlayers)
resolution : 0.1666667, 0.1666667 (x, y)
extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names : alt.1, alt.2, bio1, bio5, bio6, bio12, bio16, bio17
min values : -353, -353, -269, -59, -547, 0, 0, 0
max values : 6241, 6241, 314, 489, 258, 9916, 5043, 2159
> nlayers(stackrasters) #verifica número de camadas
[1] 8
Podemos remover a segunda camada com a função dropLayer:
> nlayers(stackrasters)
[1] 7
As camadas sobrepostas (stackraster) são globais, porém como queremos modelar ocorrências apenas na América do Sul, podemos recortar as camadas usando nosso raster ext:
> southamer<-crop(stackrasters,ext)
> southamer@data@names #verifica os nomes das camadas
[1] "alt.1" "bio1" "bio5" "bio6" "bio12" "bio16" "bio17"
Finalmente, podemos rodar o modelo com a função bioclim, que usa o algoritmo de mesmo nome para a modelagem. A função predict cria o raster para ser plotados, e a função reclassify apenas suaviza as cores no mapa plotado (ela não é necessária para o produto final).
> Mauritia_SDM<-bioclim(southamer,Mauritia)
> Mauritia_SDM<-predict(Mauritia_SDM, southamer)
> Mauritia_SDM<-reclassify(Mauritia_SDM,c(0,0.05,0.25, 0.051,0.1,0.5, 0.11,0.4,1))
> Mauritiella_SDM<-bioclim(southamer,Mauritiella)
> Mauritiella_SDM<-predict(Mauritiella_SDM, southamer)
> Mauritiella_SDM<-reclassify(Mauritiella_SDM,c(0,0.05,0.25, 0.051,0.1,0.5, 0.11,0.4,1))
E agora, os plots:
> par(mfrow=c(1,2))
> plot(Mauritia_SDM,main="Mauritia"); points(Mauritia, pch=19, cex=.1)
> plot(Mauritiella_SDM,main="Mauritiella"); points(Mauritiella, pch=19, cex=.1)
Podemos também combinar os dois modelos num mapa só"
> overlap<-Mauritia_SDM+Mauritiella_SDM
> plot(overlap/2)
> grid(); scalebar(1000, type='bar', divs=2,below="km")
