Preparando pacotes

Caso essa seja a sua primeira vez utilizando esses pacotes, talvez seja necessário instalar alguns deles:

  #Pacotes necessários
  package_list<-c("rgbif", "maptools")
 
  #Carrega ou instala os pacotes

for(i in package_list){
  if(!require(i, character.only = T)){
    install.packages(i, character=T) #realmente demora!
  }else{
    library(i, character.only = T)
  }
}
## Loading required package: rgbif
## Loading required package: maptools
## Loading required package: sp
## Checking rgeos availability: TRUE

1. Buscando e carregando dados

Primeiramente, vamos setar o diretório onde iremos trabalhar

#setando diretório
setwd("/media/wilson/MIDGARD/Documents/AARE - MNE/Aula 4 - Load and Data Cleaning/")

Neste exerício vamos importar dados diretamente do Global Biodiversity Information Facility (gbif.org) através do pacote rgbif

species <- "Euglossa annectans"
occs <- occ_search(scientificName = species, 
                   return = "data")
## Warning: `return` param in `occ_search` function is defunct as of rgbif v3.0.0, and is ignored
## See `?rgbif` for more information.
## This warning will be thrown once per R session.
occs
## Records found [126] 
## Records returned [126] 
## No. unique hierarchies [1] 
## No. media records [126] 
## No. facets [0] 
## Args [limit=500, offset=0, scientificName=Euglossa annectans, fields=all] 
## # A tibble: 126 x 103
##    key   scientificName decimalLatitude decimalLongitude issues datasetKey
##    <chr> <chr>                    <dbl>            <dbl> <chr>  <chr>     
##  1 3290… Euglossa anne…           -23.5            -45.8 gass8… 86aba78e-…
##  2 3290… Euglossa anne…           -23.8            -45.4 gass8… 86aba78e-…
##  3 3290… Euglossa anne…           -23.5            -46.2 gass8… 86aba78e-…
##  4 3290… Euglossa anne…           -23.5            -45.8 gass8… 86aba78e-…
##  5 3290… Euglossa anne…           -23.5            -45.8 gass8… 86aba78e-…
##  6 3290… Euglossa anne…           -23.5            -46.2 gass8… 86aba78e-…
##  7 3290… Euglossa anne…           -23.5            -46.2 gass8… 86aba78e-…
##  8 3290… Euglossa anne…           -23.5            -46.2 gass8… 86aba78e-…
##  9 3290… Euglossa anne…           -23.5            -46.2 gass8… 86aba78e-…
## 10 3290… Euglossa anne…           -23.5            -46.2 gass8… 86aba78e-…
## # … with 116 more rows, and 97 more variables: publishingOrgKey <chr>,
## #   installationKey <chr>, publishingCountry <chr>, protocol <chr>,
## #   lastCrawled <chr>, lastParsed <chr>, crawlId <int>,
## #   hostingOrganizationKey <chr>, basisOfRecord <chr>, occurrenceStatus <chr>,
## #   sex <chr>, taxonKey <int>, kingdomKey <int>, phylumKey <int>,
## #   classKey <int>, orderKey <int>, familyKey <int>, genusKey <int>,
## #   speciesKey <int>, acceptedTaxonKey <int>, acceptedScientificName <chr>,
## #   kingdom <chr>, phylum <chr>, order <chr>, family <chr>, genus <chr>,
## #   species <chr>, genericName <chr>, specificEpithet <chr>, taxonRank <chr>,
## #   taxonomicStatus <chr>, continent <chr>, stateProvince <chr>, year <int>,
## #   month <int>, day <int>, eventDate <chr>, modified <chr>,
## #   lastInterpreted <chr>, license <chr>, identifiers <chr>, facts <chr>,
## #   relations <chr>, isInCluster <lgl>, geodeticDatum <chr>, class <chr>,
## #   countryCode <chr>, recordedByIDs <chr>, identifiedByIDs <chr>,
## #   country <chr>, rightsHolder <chr>, identifier <chr>, locality <chr>,
## #   municipality <chr>, datasetName <chr>, gbifID <chr>, language <chr>,
## #   collectionCode <chr>, occurrenceID <chr>, type <chr>, catalogNumber <chr>,
## #   recordedBy <chr>, institutionCode <chr>, ownerInstitutionCode <chr>,
## #   name <chr>, gadm <chr>, lifeStage <chr>, http...unknown.org.rights <chr>,
## #   http...unknown.org.language <chr>, http...unknown.org.rightsHolder <chr>,
## #   fieldNumber <chr>, occurrenceRemarks <chr>, identifiedBy <chr>,
## #   higherGeography <chr>, institutionID <chr>, county <chr>,
## #   preparations <chr>, rights <chr>, coordinateUncertaintyInMeters <dbl>,
## #   institutionKey <chr>, collectionKey <chr>, georeferencedBy <chr>,
## #   georeferenceProtocol <chr>, accessRights <chr>, dateIdentified <chr>,
## #   typeStatus <chr>, X.99d66b6c.9087.452f.a9d4.f15f2c2d0e7e. <chr>,
## #   individualCount <int>, verbatimEventDate <chr>, endDayOfYear <chr>,
## #   startDayOfYear <chr>, higherClassification <chr>, verbatimLocality <chr>,
## #   waterBody <chr>, taxonID <chr>, locationID <chr>, samplingProtocol <chr>

Note que o objeto “occs” é um tipo complexo, constituído por cinco listas ou slots. Uma delas (“data”) é uma tabela especial chamada tibble (Tyddiverse) que estoca dados e metadados e pode também ser coercida à dataframe.

Checando o tamanho do dataset

Uma função básica para checar tabelas de dados é a função dim(). Ela retorna número de linhas e colunas em uma tabela / dataset.

dim(occs$data)
## [1] 126 103

Também podemos checar a lista de colunas no documento:

head(colnames(occs$data)) #Padrão Darwin Core
## [1] "key"              "scientificName"   "decimalLatitude"  "decimalLongitude"
## [5] "issues"           "datasetKey"
#Note que utilizei aqui a função head() que retorna as primeiras 10 linhas da tabela. Caso queira ver todas, apenas remova o bloco colnames(occs$data) de dentro da função head().

Salvando dados originais

Documentar todos os passos é fundamental. Neste ponto, salvar os dados originais é mandatório.

As próximas linhas de código, primeiro criam a pasta “data” (caso ela não exista) e em seguida salvam os dados ‘occs$data’ como arquivo csv:

dir.create("data")
## Warning in dir.create("data"): 'data' already exists
write.csv(occs$data, 
          "data/raw_data.csv", 
          row.names = FALSE)

2. Checando dados da espécie

Uma vez que temos os dados brutos, um passo importante é checar a consistência do nome específico.

Primeiro vamos checar quantos nomes científicos foram importados. Note que a função unique() retorna todos os nomes únicos para a coluna indicada:

sort(unique(occs$data$scientificName))
## [1] "Euglossa annectans Dressler, 1982"

Neste caso, temos apenas um nome, então estamos bem para prosseguir.

Também podemos checar outros campos do banco de dados:

sort(unique(occs$data$species))
## [1] "Euglossa annectans"
sort(unique(occs$data$speciesKey))
## [1] 9101860

Veja que para os dois comandos, foram retornados dados únicos. Assim, podemos seguir sem precisar consolidar esses dados.

Embora simples esse passos são fundamentais e pode interferir em toda modelagem. Gaste o tempo que for necessário aqui para confirmar se a espécie está corretamente anotada.

Dica importante: caso note que é necessário alterar algo em relação ao nome da espécie, nunca faça isso no objeto original, crie um objeto novo, algo como:

dados<-occs$data. Em seguida edite faça as alterações sobre o objeto dados.

3. Checando coordenadas

Uma vez que definimos que não há problemas com a espécie, podemos seguir adiante.

Vamos primeiro plotar a latitude e longitude para observar seu padrão:

plot(decimalLatitude ~ decimalLongitude, data = occs$data)

#map(, , , add = TRUE)

Plot lat ~ long

Observe como temos várias ocorrências sobrepostas entre -20S e -35W. Essas sobreposições são certamente duplicatas e precisaremos remove-las. Do mesmo modo, temos outlyers próximo à 0 graus de latitude e longitude, além de um dado entre -50S e -25W. Certamente esses dados são erros e precisaremos corrigí-los.

Primero checaremos a unicidade dos dados…

Para esse exercício não precisaremos de redundância, portanto trabalharemos apenas com registros únicos, que não se repetem tampouco em Lat ou Long.

Aplicaremos a função “levels” que retorna os níveis (elementos únicos) dos fatores (variável categórica). A função “as.factor” transforma em fator (categorias) o parâmetro que é passado, neste caso, os dados da coluna decimalLatitude.

levels(as.factor(occs$data$decimalLatitude))
##  [1] "-55"        "-23.8311"   "-23.6039"   "-23.5475"   "-23.5322"  
##  [6] "-23.5228"   "-23.423056" "-23.1864"   "-22.95831"  "-22.5425"  
## [11] "-19.90174"  "0"

Fazemos o mesmo para longitude

levels(as.factor(occs$data$decimalLongitude))
##  [1] "-56.487472" "-52.327778" "-46.9192"   "-46.8842"   "-46.6361"  
##  [6] "-46.1883"   "-45.8464"   "-45.4264"   "-43.9642"   "-43.27721" 
## [11] "-26"        "0"

Aqui fica claro que temos alguns dados imprecisos (-55 em Long ou -26 em Lat) e alguns que foram erroneamente setados como 0.

Note ainda que, embora temos checado os dados, ainda não os alteramos de forma direta, faremos isso a partir de agora

Corrigindo os dados

Primeiramente, removeremos os valores tidos como não disponíveis, ou NA.

occs.coord <- occs$data[!is.na(occs$data$decimalLatitude) 
                   & !is.na(occs$data$decimalLongitude),]

head(occs.coord)
## # A tibble: 6 x 103
##   key   scientificName decimalLatitude decimalLongitude issues datasetKey
##   <chr> <chr>                    <dbl>            <dbl> <chr>  <chr>     
## 1 3290… Euglossa anne…           -23.5            -45.8 gass8… 86aba78e-…
## 2 3290… Euglossa anne…           -23.8            -45.4 gass8… 86aba78e-…
## 3 3290… Euglossa anne…           -23.5            -46.2 gass8… 86aba78e-…
## 4 3290… Euglossa anne…           -23.5            -45.8 gass8… 86aba78e-…
## 5 3290… Euglossa anne…           -23.5            -45.8 gass8… 86aba78e-…
## 6 3290… Euglossa anne…           -23.5            -46.2 gass8… 86aba78e-…
## # … with 97 more variables: publishingOrgKey <chr>, installationKey <chr>,
## #   publishingCountry <chr>, protocol <chr>, lastCrawled <chr>,
## #   lastParsed <chr>, crawlId <int>, hostingOrganizationKey <chr>,
## #   basisOfRecord <chr>, occurrenceStatus <chr>, sex <chr>, taxonKey <int>,
## #   kingdomKey <int>, phylumKey <int>, classKey <int>, orderKey <int>,
## #   familyKey <int>, genusKey <int>, speciesKey <int>, acceptedTaxonKey <int>,
## #   acceptedScientificName <chr>, kingdom <chr>, phylum <chr>, order <chr>,
## #   family <chr>, genus <chr>, species <chr>, genericName <chr>,
## #   specificEpithet <chr>, taxonRank <chr>, taxonomicStatus <chr>,
## #   continent <chr>, stateProvince <chr>, year <int>, month <int>, day <int>,
## #   eventDate <chr>, modified <chr>, lastInterpreted <chr>, license <chr>,
## #   identifiers <chr>, facts <chr>, relations <chr>, isInCluster <lgl>,
## #   geodeticDatum <chr>, class <chr>, countryCode <chr>, recordedByIDs <chr>,
## #   identifiedByIDs <chr>, country <chr>, rightsHolder <chr>, identifier <chr>,
## #   locality <chr>, municipality <chr>, datasetName <chr>, gbifID <chr>,
## #   language <chr>, collectionCode <chr>, occurrenceID <chr>, type <chr>,
## #   catalogNumber <chr>, recordedBy <chr>, institutionCode <chr>,
## #   ownerInstitutionCode <chr>, name <chr>, gadm <chr>, lifeStage <chr>,
## #   http...unknown.org.rights <chr>, http...unknown.org.language <chr>,
## #   http...unknown.org.rightsHolder <chr>, fieldNumber <chr>,
## #   occurrenceRemarks <chr>, identifiedBy <chr>, higherGeography <chr>,
## #   institutionID <chr>, county <chr>, preparations <chr>, rights <chr>,
## #   coordinateUncertaintyInMeters <dbl>, institutionKey <chr>,
## #   collectionKey <chr>, georeferencedBy <chr>, georeferenceProtocol <chr>,
## #   accessRights <chr>, dateIdentified <chr>, typeStatus <chr>,
## #   X.99d66b6c.9087.452f.a9d4.f15f2c2d0e7e. <chr>, individualCount <int>,
## #   verbatimEventDate <chr>, endDayOfYear <chr>, startDayOfYear <chr>,
## #   higherClassification <chr>, verbatimLocality <chr>, waterBody <chr>,
## #   taxonID <chr>, locationID <chr>, samplingProtocol <chr>
print(paste("com NA:", nrow(occs$data))) #com NA
## [1] "com NA: 126"
print(paste("sem NA:", nrow(occs.coord))) #sem NA
## [1] "sem NA: 107"

Após remover os “NA” podemos salvar os dados em um novo objeto (chamado aqui “gbif”) que será um dataframe com três colunas: species, decimalLongitude e decimalLatitude. Note que a ordem desses campos é importate. Devem corresponder à nome da espécie, valor para o eixo X e valor para o eixo Y + PS: alterei a ordem dos campos em relação ao vídeo, agora estão como deveriam estar

Note que aqui não utilizaremos os dados adicionais da tabela, por isso nos damos ao luxo de utilizar apenas as três colunas.

#Criando uma tabela nova apenas com spp, long e lat.
gbif<-data.frame("species" = occs.coord$species,
                 "decimalLongitude" = occs.coord$decimalLongitude,
                  "decimalLatitude" = occs.coord$decimalLatitude
                 )
head(gbif) #primeiras 10 linhas
##              species decimalLongitude decimalLatitude
## 1 Euglossa annectans         -45.8464        -23.5322
## 2 Euglossa annectans         -45.4264        -23.8311
## 3 Euglossa annectans         -46.1883        -23.5228
## 4 Euglossa annectans         -45.8464        -23.5322
## 5 Euglossa annectans         -45.8464        -23.5322
## 6 Euglossa annectans         -46.1883        -23.5228
tail(gbif) #últimas 10 linhas
##                species decimalLongitude decimalLatitude
## 102 Euglossa annectans        -56.48747       -23.42306
## 103 Euglossa annectans        -52.32778       -22.54250
## 104 Euglossa annectans        -26.00000       -55.00000
## 105 Euglossa annectans        -52.32778       -22.54250
## 106 Euglossa annectans        -52.32778       -22.54250
## 107 Euglossa annectans          0.00000         0.00000

Há formas alternativas de visualizar os dados:

#vizualizando os dados no Rstudio
View(gbif)

Podemos contar os registros:

#contando quantas linhas há na tabela gbif
nrow(gbif)
## [1] 107

E quantos restaram após remover os dados não disponíveis:

#Numero de linhas no dado original
nrow(occs$data)
## [1] 126
#Original menos o modificado
nrow(occs$data)-nrow(gbif)
## [1] 19

Podemos novamente visualizar as latitudes.

#Visual checking
plot(gbif$decimalLongitude)

plot(gbif$decimalLatitude)

[Plotes de Latitude e longitude. Note que nestes plotes há uma diferença entre plotar Latitude contra longitude. Aqui podemos ver a repetição dos valores (em Y) ao longo dos registros (em X).]

Checando os dados no mapa:

As próximas linhas de código vão importar os dados de vetor do maptools de um conjunto chamado wrld_simpl que possue as fronteiras dos países do mundo.

#data from maptools
data("wrld_simpl")

As próximas linhas efetivamente plotam o mapa.

#plotando o mapa
plot(wrld_simpl, xlim=c(-90, 90), ylim=c(-60,60), axes=T, col='olivedrab3',bg='lightblue')
#Para aproximar o plot mude xlim para: xlim=c(-60, -55)
#                           ylim para: ylim=c(-40,-15)
points(gbif[ ,c(2,3)], col="#FF0000")

#Note que em pontos, as colunas 2, 3 devem representar Long e Lat, nessa ordem. Vídeo possui ordem desses fatores trocada.

Distribuição dos pontos de coleta sobre o sul do Brasil.

Note que é possível mudar o intervalo de plot alterando os parâmetros xlim e ylim.

Remover os pontos com valor zero ou fora do “range”

Note que as linhas a seguir executam uma série de códigos para manter os pontos entre -40 e -18 de Lat(Sul e Norte) e -40 e -60 Long (Leste - Oeste)

#Mantendo dados diferentes de 0 em lat e long
gbif <- 
gbif[gbif$decimalLatitude != 0 & gbif$decimalLongitude != 0, ]
#Mantendo registros em Lat maior que -40S **ou** maior q long -40W
gbif <- 
gbif[gbif$decimalLatitude > -40 | gbif$decimalLongitude > -40, ]
#Mantendo dados menores que -18S **ou** -60W
gbif <- 
gbif[gbif$decimalLatitude < -18 | gbif$decimalLongitude < -60, ]
#Mantendo valores diferentes de valores específicos 
gbif <- 
gbif[gbif$decimalLatitude != 55.000 -18 & gbif$decimalLongitude != -26.00000, ]

Podemos agora, plotar os pontos apenas para ver sua distribuição:

#checking again
plot(gbif$decimalLongitude)

plot(gbif$decimalLatitude)

Plots de latitude e longitude.

Note que plotamos latitude ou longitude em Y e seus registros ou linhas da tabela em X. Logo cada ponto em x corresponde à uma linha e cada valor correspondente em y uma latitude ou longitude. Note como os dados se repetem. É necessário remover as duplicatas.

Removendo duplicatas

Remover duplicatas no R é fácil. Basta usar a função nativa “unique”. Essa função retarna apenas valores únicos. Lembre-se de executar essa função em um objeto novo para não sobrescrever nada.

gbif.unique <- unique(gbif)
gbif.unique
##                species decimalLongitude decimalLatitude
## 1   Euglossa annectans        -45.84640       -23.53220
## 2   Euglossa annectans        -45.42640       -23.83110
## 3   Euglossa annectans        -46.18830       -23.52280
## 36  Euglossa annectans        -46.63610       -23.54750
## 89  Euglossa annectans        -43.96420       -19.90174
## 91  Euglossa annectans        -46.88420       -23.18640
## 92  Euglossa annectans        -46.91920       -23.60390
## 98  Euglossa annectans        -43.27721       -22.95831
## 100 Euglossa annectans        -56.48747       -23.42306
## 103 Euglossa annectans        -52.32778       -22.54250

Plotando no mapa novamente

plot(wrld_simpl, xlim=c(-60, -55), ylim=c(-40,-15), axes=T, col='olivedrab3',bg='lightblue')
points(gbif.unique[ ,c(2,3)], col="#FF0000")#lembre-se, pontos devem ser long e lat, nessa ordem. Ou seja, as colunas em gbif.unique devem ser "species", "decimalLong", "decimalLat".

#Salvando os dados trabalhados

Uma vez que os dados foram corrigidos, podemos salvá-los como entrada para os algoritmos de modelagem.

#Creating an input for modelling
input <- gbif.unique

#saving the input table
write.table(input, "input.coordinates.txt", quote=F, row.names=F, sep = "\t")

#Verificando se os dados foram salvos:

Vamos agora importar os dados para checar se foram salvos corretamente:

input.test <- read.csv("input.coordinates.txt", header = T, sep = "\t")
head(input.test)
##              species decimalLongitude decimalLatitude
## 1 Euglossa annectans         -45.8464       -23.53220
## 2 Euglossa annectans         -45.4264       -23.83110
## 3 Euglossa annectans         -46.1883       -23.52280
## 4 Euglossa annectans         -46.6361       -23.54750
## 5 Euglossa annectans         -43.9642       -19.90174
## 6 Euglossa annectans         -46.8842       -23.18640
nrow(input.test)
## [1] 10
plot(input.test[,2:3])

[Se tudo correu bem, o gráfico deve ser como o da figura, mostrando como as ocorrências se distribuem em lat e long]

Como ultima checagem, podemos olhar as ocorrências importadas em um mapa:

plot(wrld_simpl, xlim=c(-60, -55), ylim=c(-40,-15), axes=T, col='olivedrab3',bg='lightblue')
points(input.test[ ,c(2,3)], col="#FF0000")#lembre-se, pontos devem ser long e lat, nessa ordem. Ou seja, as colunas em gbif.unique devem ser "species", "decimalLong", "decimalLat".

Pontos plotados sobre o mapa cobrindo a porção sul da Mata Altântica.

Espero que tenha aproveitado o tutorial. Vimos um pouco sobre como os erros dos bancos de dados podem “sujar” nossas análises, e que esses podem estar incompletos. Em breve veremos como aglutinar dados de diferentes bancos de dados para as análises.

Bons estudos!

Prof. Dr. Wilson Frantine-Silva
Pesquisador PNPD - ORCID
Universidade Estadual do Norte Fluminense Darcy Ribeiro - LCA.
Programa de Pós-Graduação em Ecologia e Recursos Naturais