Aplicações

7.0.1 Pontos de incêndio nas florestas brasileiras

Para a segunda aplicação, será utilizada uma base de dados do kaggle com os pontos e a quantidade de incêndios nas florestas brasileiras. Para isso, serão instalados os pacotes abaixo.

options(repos = list(CRAN="http://cran.rstudio.com/"))
options("install.lock"=FALSE)

install.packages(c('sf','mapview'))

library(sf)
library(mapview)
library(ggplot2)     # graphic
library(tidyverse)
library(dplyr)       # data manipulation
library(readr)       # read files

O dataset foi retirado do Kaggle e o link está apresentado no code chunk abaixo.

# Data frame com os incêndios e os estados brasileiros
# fonte: https://www.kaggle.com/code/tianjinglei/starter-forest-fires-in-brazil-41eec699-a/input

forest_fires <- read.csv("Data/amazon.csv")
head(forest_fires)

##   year state   month number       date
## 1 1998  Acre Janeiro      0 1998-01-01
## 2 1999  Acre Janeiro      0 1999-01-01
## 3 2000  Acre Janeiro      0 2000-01-01
## 4 2001  Acre Janeiro      0 2001-01-01
## 5 2002  Acre Janeiro      0 2002-01-01
## 6 2003  Acre Janeiro     10 2003-01-01

View(forest_fires)
str(forest_fires)

## 'data.frame':    6454 obs. of  5 variables:
##  $ year  : int  1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 ...
##  $ state : chr  "Acre" "Acre" "Acre" "Acre" ...
##  $ month : chr  "Janeiro" "Janeiro" "Janeiro" "Janeiro" ...
##  $ number: num  0 0 0 0 0 10 0 12 4 0 ...
##  $ date  : chr  "1998-01-01" "1999-01-01" "2000-01-01" "2001-01-01" ...

Analisando a base de dados, é possível perceber que houveram problemas na leitura de estados que possuem acento no nome, como Pará e Piauí. Para resolver isso:

# Problema com a palavra "Pará" porque tem um acento e o R não está conseguindo ler

Encoding(forest_fires$state) <- "latin1"
forest_fires$state[forest_fires$state == 'Par�'] <- 'Para' #Ajuste

# Piauí tambémdeu problema
forest_fires$state[forest_fires$state == 'Piau'] <- 'Piaui' #Ajuste
View(forest_fires)

Para analisar a quantidade de incêndios ao longo dos anos nos estados brasileiros, será realizada a soma deles em cada estado:

# Agrupando a soma dos incêndios por estado

total_fires_by_state <- forest_fires %>%
  group_by(state) %>%
  summarise(sum = sum(number))

# Visualizando a base de dados criada com a soma por estado
View(total_fires_by_state)

# Sumarização dos dados
summary(total_fires_by_state)

##     state                sum       
##  Length:23          Min.   : 3237  
##  Class :character   1st Qu.:21059  
##  Mode  :character   Median :25129  
##                     Mean   :30388  
##                     3rd Qu.:37750  
##                     Max.   :96246

Para visualizar a quantidade de incêndios por estado, será realizado um gráfico de barras:

# Gráfico de barras com a soma de incêncios por estado
ggplot(total_fires_by_state, aes(x = state, y = sum)) +
  geom_bar(stat='identity') +
  coord_flip()

Para poder realizar a criação de mapas, foi criado um arquivo com as coordenadas (retiradas do Google Maps) de cada um dos estados brasileiros em WGS84:

# Dataframe com as coordenadas lat e long em WGS84 (fonte: google maps)
coord <- read.table("Data/coordenadas.txt", sep=",",dec=".",header = TRUE)
View(coord)

Para poder analisar a quantidade de incêndios no mapa, será realizada a junção entre a tabela com a quantidade de incêndio por estado e a tabela de coordenadas:

# Juntando os dataframes
joindata <- inner_join(total_fires_by_state,coord,by="state")
View(joindata)

Para que as coordenadas da tabela sejam tratadas como simple features, será utilizada a função st_as_sf() para criar os pontos das coordenadas de acordo com a latitude e longitude:

# Junção entre a tabela com as coordenadas e a quantidade de incêndios:
joindata2 <- inner_join(total_fires_by_state,coord,by="state") %>%
  st_as_sf(coords=c('long','lat'), crs=4326) # criação dos pontos com as coordenadas e inserção do CRS de referência

Para poder visualizar cada um dos pontos com os limites geográficos dos estados brasileiros, será importado um shapefile do IBGE. Sabe-se que os shapefiles do IBGE como CRS o SIRGAS 2000 e isso será confirmado a partir da função ´st_crs()´.

# Importando shapefile (fonte: ibge)

shape <- read_sf('Data/Shapefiles/BR_UF_2020.shp') # por ser do ibge, utiliza crs = SIRGAS2000 (EPSG code: 4674), precisamos mudar para WGS84
st_crs(shape)

## Coordinate Reference System:
##   User input: SIRGAS 2000 
##   wkt:
## GEOGCRS["SIRGAS 2000",
##     DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
##         ELLIPSOID["GRS 1980",6378137,298.257222101,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore."],
##         BBOX[-59.87,-122.19,32.72,-25.28]],
##     ID["EPSG",4674]]

View(shape)

Como o objetivo é a visualização da ´joindata´ e do shapefile em um mesmo mapa, ambos precisam estar com o mesmo CRS. Dessa forma, o shapefile do IBGE será transformado em WGS84 a partir da função `st_transform:

shape2 <- st_transform(shape, 4326) # Mudando crs para WGS84 (EPSG code: 4326)
st_crs(shape2)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

Antes de criarmos os mapas, iremos checr se os CRS são os mesmos:

# Checando SGR
st_crs(joindata2) 

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

st_crs(shape2)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

O primeiro mapa será um mapa interativo:

# Mapa Interativo com os pontos no mapa do Brasil
mapview(shape2, col.regions = 'green') +
  mapview(joindata2)

O segundo mapa levará em conta a quantidade de incêndio em cada estado. Será criada uma escala de cor de acordo com a quantidade de incêndios. Assim, poderemos analisar quais estados apresentam maior quantidade de incêndios nas florestas ao longo dos anos.

# Mapa com cor de ponto diferente de acordo com a quantidade de incêndios
ggplot() + 
  geom_sf(data = shape2) + 
  geom_sf(data = shape2, colour = "green", fill = NA) + 
  geom_point(data = joindata, mapping = aes(x = long, y = lat, colour = sum, size = 3)) +
    scale_colour_gradient2(
    low = "yellow", 
    mid = "orange", 
    high = "red",
    midpoint = 30388) +
  coord_sf() +
  ggtitle("Pontos de incêndio no Brasil")

7.0.2 Relação entre quantidade de algas e de peixes em Tampa Bay

Para a aplicação será reproduzido o exemplo do TBEP R Training. O shapefile “sgdat.shp” são dados da cobertura de algas marinhas em Tampa Bay em 2016. As “features” são as linhas do vetor e os “fields” são as colunas, ou melhor, atributos (“OBJECT ID” e “FLUCCS”). O SGR do arquivo é WGS 84. A coluna “geometry” armazena os dados espaciais (longitude e latitude).

Esse é o passo a passo de como importar um shapefile. Porém, muitas vezes não possuímos um shapefile e queremos criar um a partir de um dataframe. Para isso, é necessário que o dataframe inclua as coordenadas geográficas (longitude e latitude) e que tenhamos conhecimento do SGR. O dataframe ´fishdat´ possui as características dos peixes encontrados e o statloc apresenta a localização deles. O passo a passo será realizado abaixo.

# dados da presença de peixes em Tampa Bay
fishdat <- read.csv("Data/fishdat.csv")

#localização geográfica dos peixes
statloc <- read.csv("Data/statloc.csv")

# estrutura dos dados
str(fishdat)

## 'data.frame':    2844 obs. of  12 variables:
##  $ OBJECTID     : int  1550020 1550749 1550750 1550762 1550828 1550838 1550842 1551131 1551311 1551335 ...
##  $ Reference    : chr  "TBM1996032006" "TBM1996032004" "TBM1996032004" "TBM1996032207" ...
##  $ Sampling_Date: chr  "1996-03-20" "1996-03-20" "1996-03-20" "1996-03-22" ...
##  $ yr           : int  1996 1996 1996 1996 1996 1996 1996 1996 1996 1996 ...
##  $ Gear         : int  300 22 22 20 160 300 300 300 300 22 ...
##  $ ExDate       : chr  "2018-04-12 10:27:38" "2018-04-12 10:25:23" "2018-04-12 10:25:23" "2018-04-12 10:25:23" ...
##  $ Bluefish     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Common.Snook : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Mullets      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Pinfish      : int  0 54 0 80 0 0 0 0 1 1 ...
##  $ Red.Drum     : int  0 0 1 0 4 0 0 0 0 0 ...
##  $ Sand.Seatrout: int  1 0 0 0 0 1 5 66 0 0 ...

str(statloc)

## 'data.frame':    2173 obs. of  3 variables:
##  $ Reference: chr  "TBM1996032006" "TBM1996032004" "TBM1996032207" "TBM1996042601" ...
##  $ Latitude : num  27.9 27.9 27.9 28 27.9 ...
##  $ Longitude: num  -82.6 -82.6 -82.5 -82.7 -82.6 ...

Para isso, utilizaremos a função st_as_sf() para transformar o dataframe em um objeto sf. Primeiramente, precisamos juntar os dois datasets (fishdat e statloc) e dizer qual coluna que possui os dados da geometria (latitude e longitude). Além disso, é necessário dizer qual o SGR e, além disso, precisamos garantir que ambos datasets possuam o mesmo SGR. Por enquanto, podemos fazer um “chute calibrado” que é o WGS84.

#juntando os dois dataframes
library(dplyr)
joindata <- left_join(fishdat,statloc,by="Reference")

#criando o objeto de dados espaciais
joindata <- st_as_sf(joindata, coords=c('Longitude','Latitude'), crs = st_crs(sgdat))

#tipo de objeto sf
str(joindata)

## Classes 'sf' and 'data.frame':   2844 obs. of  13 variables:
##  $ OBJECTID     : int  1550020 1550749 1550750 1550762 1550828 1550838 1550842 1551131 1551311 1551335 ...
##  $ Reference    : chr  "TBM1996032006" "TBM1996032004" "TBM1996032004" "TBM1996032207" ...
##  $ Sampling_Date: chr  "1996-03-20" "1996-03-20" "1996-03-20" "1996-03-22" ...
##  $ yr           : int  1996 1996 1996 1996 1996 1996 1996 1996 1996 1996 ...
##  $ Gear         : int  300 22 22 20 160 300 300 300 300 22 ...
##  $ ExDate       : chr  "2018-04-12 10:27:38" "2018-04-12 10:25:23" "2018-04-12 10:25:23" "2018-04-12 10:25:23" ...
##  $ Bluefish     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Common.Snook : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Mullets      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Pinfish      : int  0 54 0 80 0 0 0 0 1 1 ...
##  $ Red.Drum     : int  0 0 1 0 4 0 0 0 0 0 ...
##  $ Sand.Seatrout: int  1 0 0 0 0 1 5 66 0 0 ...
##  $ geometry     :sfc_POINT of length 2844; first list element:  'XY' num  -82.6 27.9
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:12] "OBJECTID" "Reference" "Sampling_Date" "yr" ...

#checando SGR
st_crs(joindata)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

st_crs(sgdat) 

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

Caso seja necessário modificar a projeção, utiliza-se a função `st_transform(). Nesse caso, não precisamos modificar já que o shapefile (sgdat) tem o mesmo SGR do que estamo querendo criar.

Agora, iniciaremos a análise geoespacial dos dados. Inicialmente, iremos dar uma olhada geral para entender qual os dados que estamos lidando.

O padrão é que a função ´plot()´ plote todas as feições. Para plotar somente a geometria, utiliza-se st_geometry().

plot(st_geometry(joindata)) 

plot(joindata)

plot(sgdat)

Conforme observamos o shapefile “sgdat” com os dados das algas marinhas e o “joindata” com os dados do posicionamento de peixes, é possível verificar que existem áreas de intersecção entre ambos. Para analisar novamente, iremos plotar somente a geometria de ambos:

plot(joindata$geometry)

plot(sgdat$geometry)

Vamos filtrar somente os dados dos peixes do ano de 2016:

filt_data <- joindata %>%
  filter(yr == 2016)
plot(st_geometry(filt_data))

Agora, verificaremos quantos peixes foram vistos nos mesmos locais em que encontraram-se algas marinhas em 2016. Ou seja, iremos selecionar as localizações que possuem ambos dados. Para isso, iremos utilizar o código abaixo:

fish_crop <- filt_data[sgdat, ]

plot(fish_crop$geometry)

O que foi realizado até agora é somente a intersecção da geometria de ambos datasets. Portanto, agora realizaremos a intersecção de ambos dados, incluindo atributos:

fish_int <- st_intersection(filt_data, sgdat)

plot(st_geometry(fish_int))

View(fish_int)

É possível utilizar ferramentas do tidyverse. Abaixo, iremos fazer a soma de todos os Pinfish foram pegos em 2016:

fish_cnt <- fish_int %>% 
  group_by(FLUCCS) %>% 
  summarise(
    cnt = sum(Pinfish)
  ) 

fish_cnt

## Simple feature collection with 2 features and 2 fields
## Geometry type: MULTIPOINT
## Dimension:     XY
## Bounding box:  xmin: -82.7182 ymin: 27.82623 xmax: -82.53237 ymax: 28.02418
## Geodetic CRS:  WGS 84
## # A tibble: 2 × 3
##   FLUCCS   cnt                                                          geometry
##   <chr>  <int>                                                  <MULTIPOINT [°]>
## 1 9113    1559 ((-82.53352 27.92907), (-82.53617 27.91043), (-82.53535 27.91212…
## 2 9116    4766 ((-82.55893 27.96612), (-82.5588 27.96633), (-82.55797 27.96632)…

Além de realizar a soma nos atributos numéricos (quantidade de Pinfish), também é realizada nos atributos geométricos (latitude e longitude). Conforme apresentado na tabela anterior, existe uma maior quantidade de Pinfishs em áreas onde existe maior quantidade de algas marinhas (FLUCCS=9116). É possível realizar um gráfico em relação às duas categorias de cobertura de algas marinhas (´9113´: desigual, ´9116´: contínua).

ggplot(fish_cnt, aes(x = FLUCCS, y = cnt)) + 
  geom_bar(stat = 'identity', fill='navyblue') 

Agora será realizada a confecção de mapas. Utilizaremos os pacotes ggplot2 inicialmente:

ggplot() + 
  geom_sf(data = sgdat, fill = 'green') + 
  geom_sf(data = joindata) 

Agora, para criar um mapa interativo para selecionar e dar zoom nos dados, utilizaremos o pacote mapview:

mapview(sgdat, col.regions = 'green') +
  mapview(joindata, zcol = 'Gear')