Tutorial - como gerar os gráficos boxplot no R a partir do RSÁgua - Fepam

Primeiro acesso

Olá! Nesse tutorial você irá aprender a replicar tudo que será exibido nessa página!

Para isso, você precisa ter instalado em seu computador 2 programas: o R e o RStudio.
O R é a linguagem de programação e o RStudio a interface gráfica,

Você pode baixar o R e o RStudio através do link https://posit.co/download/rstudio-desktop/

1. Se é o seu primeiro acesso, clique na caixinha Code no canto superior direito para baixar o código no formato .Rmd, ou seja, RMarkdown.
2. Abra o RStudio
3. Clique em File -> Open file -> e selecione esse arquivo .Rmd que você baixou.
4. Você deve chegar a uma tela semelhante a da figura abaixo

Na primeira vez em que for rodar os códigos aqui no R, será solicitado que instale pacotes. Quando aparecerem caixas de diálogo, clique em Yes.

Pode demorar um pouco esse processo. Para ter uma melhor visualização desse tutorial pressione as teclas Ctrl+Shift+K ou clique no botão Knit acima. Recomendo manter uma aba do navegador com o arquivo tutorial_R.html aberto.

Já podemos rodar o código pela primeira vez!
Para isso, aperte Ctrl+Alt+R ou clique em Run no menu superior e depois em Run All:

Enquanto os pacotes estão sendo instalados, prossiga com a leitura. Deixei alguns conteúdos nas Referências que irão complementar esse documento.

Quando quiser renderizar o código para ter a mesma visualização como a do formato .html use o comando de Knit, o Ctrl+Shift+K, ou caso esteja com um gráfico aberto, clique em Viewer na janela inferior direita.

1 Comandos mais utilizados

Ctrl+Alt+R -> executa TODOS os códigos de uma vez só.

Knit ou Ctrl+Shift+K -> roda TODOS os códigos de uma vez só e renderiza a visualização em formatos .html ou .pdf.

Ctrl+Enter -> executa um “parágrafo” de códigos.

Alt+Enter -> executa o “parágrafo” de códigos mas sem mover o cursor.

Ctrl+Shift+C -> para comentar ou descomentar um intervalo de linhas.

As linhas de código que começam com # são comentários. Elas não impactam o resultado final do código.

Ctrl+Alt+I -> cria um novo code chunk, isto é, “pedaço de código”, é a área executável das funções que iremos utilizar.

Ctrl+I -> Indenta o código (ajusta a tabulação entre as linhas, facilita a leitura do código).

A lista completa de atalhos pode ser conferida no menu superior Tools -> Keyboard Shortcuts Help, ou através do atalho Alt+Shift+K.

2 Passo-a-passo

Importante: Com esse arquivo .Rmd você deve chegar aos mesmos resultados que presentes no arquivo tutorial_R.html.

Vou deixar comentado os pedaços de código que necessitam de alteração obrigatória com a nomenclatura alteracao_ seguido de um número.

Para encontrar facilmente os locais que precisam ser ajustados, basta utilizar o comando Ctrl+F e procurar por alteracao_ e o respectivo número listado abaixo.

1. Renomear a planilha a ser importada pra facilitar a leitura.
   • nesse caso renomeei a plan gerada pelo RSÁgua de Dados_20221222143305.xls para plan_litoral_medio.xls.
   • você pode usar Ctrl+F e substituir todas as ocorrências de plan_litoral_medio para o nome que você deu à planilha.
2. Ajustar o caminho do arquivo. (alteracao_1)
   • nesse exemplo estou importando diretamente da pasta de downloads do meu computador: "C:/Users/Léo/Downloads/plan_litoral_medio.xls"
   • Futuramente quero implementar a integração com o Google Sheets.
3. Verificar o intervalo de dados da planilha. As últimas 10 linhas costumam ser dos sumários, mas o R não entende isso.
   • verificar até qual linha a coluna Índice tem valores registrados, neste caso o Índice vai até 98 (+1 linha de cabeçalho) = linha 99
   • Alterar na parte de importação o intervalo para compreender somente os dados (alteracao_2)
   • Estou tentando descobrir como economizar essa etapa.

3 Pacotes necessários

No code chunk abaixo serão instalados todos os pacotes necessários. Cada pacote contém funções específicas que auxiliam no desenvolvimento da atividade, seja a importação do dado (readr, readxl), na manipulação dos dados (limpeza/tratamento com janitor e lubridate), na visualização (ggplot2, rmarkdown, kableExtra) ou pra praticamente tudo (tidyverse).

pacman::p_load(
  # ETL (extract, transform, load)
  janitor, readr, readxl, lubridate,
  dplyr, 
  tidyverse, 
  glue,
  # Visualização
  ggplot2,
  GGally, 
  rmarkdown, 
  knitr,
  kableExtra,
  ggbeeswarm, 
  ggtext
  # bookdown
)
#googlesheets4

Para acessar a página de ajuda dos pacotes ou das funções, basta usar um ? antes do que se está buscando.

No code chunk acima nós utilizamos uma função p_load que está dentro do pacote pacman. Para isso utilizamos o comando pacman::p_load, isto é, estamos dizendo para o R que dentro do pacote pacman quero usar especificamente a função p_load.

No exemplo abaixo estou verificando os argumentos dessa função. Na sequência abro o menu de ajuda do pacote tidyverse, um dos mais importantes da comunidade R. Outra maneira possível seria utilizando help(package = 'pacman').

?help

4 Importação dos dados

Utilizando o code chunk abaixo você estará importando a mesma planilha que utilizei. Salvei esse arquivo na nuvem, via GitHub. Estou importando-a utilizando o link para acessá-la online.

plan_litoral_medio <- read_delim(
  "https://raw.githubusercontent.com/leonardofwink/tutorial_R_fepam/main/planilha_rsagua/plan_litoral_medio%20-%20Dados_Ajustados.tsv", 
  delim = "\t", 
  escape_double = FALSE,
  col_types = cols(
    ÍNDICE = col_number(),
    `CÓD. ESTAÇÃO` = col_character(),
    LATITUDE = col_number(),
    LONGITUDE = col_number(),
    `BACIA HIDROGRÁFICA` = col_character(),
    `RECURSO HÍDRICO` = col_character(),
    REGIÃO = col_character(),
    MUNICÍPIO = col_character(),
    AMBIENTE = col_character(),
    `DATA COLETA` = col_date(format = "%m/%d/%Y"),
    `HORA COLETA` = col_time(format = "%H:%M:%S"),
    `CHUVA 24H` = col_character(),
    ALCALINIDADE = col_number(),
    ALUMÍNIO = col_number(),
    CÁDMIO = col_number(),
    CHUMBO = col_number(),
    CLORETO = col_number(),
    `CLOROFILA A` = col_number(),
    COBRE = col_number(),
    `COLIFORMES TERMOTOLERANTES` = col_number(),
    `COLIFORMES TOTAIS` = col_number(),
    CONDUTIVIDADE = col_number(),
    `CROMO TOTAL` = col_number(),
    `DEMANDA BIOQUÍMICA DE OXIGÊNIO` = col_number(),
    `DEMANDA QUÍMICA DE OXIGÊNIO` = col_number(),
    `ESCHERICHIA COLI` = col_number(),
    FERRO = col_number(),
    `FITOPLANCTON  - CIANOBACTÉRIAS` = col_number(),
    `FOSFATO ORTO` = col_number(),
    `FÓSFORO TOTAL` = col_number(),
    MANGANÊS = col_number(),
    `MERCÚRIO EM MICROGRAMA POR LITRO (UG/L)` = col_number(),
    NÍQUEL = col_number(),
    NITRATO = col_number(),
    `NITROGÊNIO AMONIACAL` = col_number(),
    `NITROGÊNIO ORGÂNICO` = col_number(),
    `NITROGÊNIO TOTAL KJELDAHL` = col_number(),
    `OXIGÊNIO DISSOLVIDO` = col_number(),
    PH = col_number(),
    `PROFUNDIDADE COLETA` = col_number(),
    `PROFUNDIDADE TOTAL` = col_number(),
    SALINIDADE = col_number(),
    `SÓLIDOS DISSOLVIDOS TOTAIS` = col_number(),
    `SÓLIDOS SUSPENSOS TOTAIS` = col_number(),
    `SÓLIDOS TOTAIS` = col_number(),
    `TEMPERATURA DA ÁGUA` = col_number(),
    `TEMPERATURA DO AR` = col_number(),
    `TRANSPARÊNCIA DA ÁGUA` = col_number(),
    TURBIDEZ = col_number(),
    `VAZÃO RECURSO HÍDRICO` = col_number(),
    ZINCO = col_number()
  ),
  locale = locale(
    date_names = "pt", 
    decimal_mark = ",",
    grouping_mark = ""
  ),
  trim_ws = TRUE
) %>%
  janitor::clean_names() %>%
  slice(
    1:(n() - 8) #retirando as ultimas 8 linhas
  ) %>% 
  rename(
    e_coli = escherichia_coli,
    dbo = demanda_bioquimica_de_oxigenio,
    mercurio = mercurio_em_micrograma_por_litro_ug_l
  ) %>%
  mutate(
    municipio = str_to_title(municipio),
    data_coleta = ymd(data_coleta),
    ano_coleta = year(data_coleta),
  ) %>%
  dplyr::select( #reordenando as colunas
    c(1:10),
    ano_coleta,
    everything()
  )

Você irá fazer as alterações nesse code chunk abaixo.
Será necessário que você altere o primeiro argumento da função read_excel, que é o caminho, ou seja: dizer ao R onde que a planilha está.

plan_litoral_medio_alterar <- read_excel(
  "C:/Users/Léo/Downloads/plan_litoral_medio.xls", #alteracao_1
  sheet = "Dados_Ajustados", 
  col_types = c(
    "numeric", "text", "numeric", "numeric", "text", 
    "text", "text", "text", "text", "date", 
    "date", "text", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", 
    "numeric"),
  range = "A1:AY99", #alteracao_2
  trim_ws = TRUE
) %>% 
  janitor::clean_names() %>% 
  rename(
    e_coli = escherichia_coli,
    dbo = demanda_bioquimica_de_oxigenio,
    mercurio = mercurio_em_micrograma_por_litro_ug_l
  ) %>% 
  mutate(
    municipio = str_to_title(municipio),
    data_coleta = ymd(data_coleta),
    ano_coleta = year(data_coleta),
    # hora_coleta = parse_datetime(hora_coleta),
  ) %>%
  dplyr::select( #reordenando as colunas
    c(1:10),
    ano_coleta,
    everything()
  )

5 Sumários estatísticos

sumario <- plan_litoral_medio %>%
  dplyr::select(cod_estacao, oxigenio_dissolvido, ano_coleta) %>%
  filter(ano_coleta > "2013" &
           ano_coleta <= "2025") %>%
  group_by(cod_estacao) %>%
  summarize(
    min =
      min(oxigenio_dissolvido,
          na.rm = TRUE),
    q1 =
      quantile(oxigenio_dissolvido, 0.20,
               na.rm = TRUE),
    median =
      median(oxigenio_dissolvido,
             na.rm = TRUE),
    mean =
      mean(oxigenio_dissolvido,
           na.rm= TRUE),
    q3 =
      quantile(oxigenio_dissolvido, 0.80,
               na.rm = TRUE),
    max =
      max(oxigenio_dissolvido,
          na.rm = TRUE))

sumario

output

## # A tibble: 6 × 7
##   cod_estacao   min    q1 median  mean    q3   max
##   <chr>       <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 87332500     0     3.3    4.94  4.75  6.73  9.23
## 2 87420130     8.14  8.58   9.06  9.15  9.67 10.6 
## 3 87420150     6.7   8.19   8.41  8.49  8.64  9.83
## 4 87420350     5.93  6.96   8.34  8.07  9.23  9.64
## 5 87420500     4.65  6.43   7.25  7.50  8.73 11.0 
## 6 87510010     2.32  3.72   5.9   5.53  7.30  8.24

5.1 Visualização da estrutura de dados da planilha

Conferir se o tipo dos dados está correto, data em formato de data (date/dttm), código da estação como character (chr), valor medido dos parâmetros como double (dbl).

glimpse(plan_litoral_medio)

output

## Rows: 98
## Columns: 52
## $ indice                      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,…
## $ cod_estacao                 <chr> "87332500", "87332500", "87332500", "87332…
## $ latitude                    <dbl> -29.91447, -29.91447, -29.91447, -29.91447…
## $ longitude                   <dbl> -50.31819, -50.31819, -50.31819, -50.31819…
## $ bacia_hidrografica          <chr> "Litoral Médio", "Litoral Médio", "Litoral…
## $ recurso_hidrico             <chr> "Lagoa dos Barros", "Lagoa dos Barros", "L…
## $ regiao                      <chr> "Litoral", "Litoral", "Litoral", "Litoral"…
## $ municipio                   <chr> "Osorio", "Osorio", "Osorio", "Osorio", "O…
## $ ambiente                    <chr> "Lêntico", "Lêntico", "Lêntico", "Lêntico"…
## $ data_coleta                 <date> 2016-06-29, 2016-09-26, 2016-12-14, 2017-…
## $ ano_coleta                  <dbl> 2016, 2016, 2016, 2017, 2017, 2017, 2017, …
## $ hora_coleta                 <time> 10:50:00, 11:05:00, 10:21:00, 11:00:00, 1…
## $ chuva_24h                   <chr> "AUSENTE", "AUSENTE", "AUSENTE", "AUSENTE"…
## $ alcalinidade                <dbl> 35.0, 11.0, 7.7, 42.5, 22.1, 14.9, 19.4, 3…
## $ aluminio                    <dbl> NA, NA, NA, NA, NA, NA, NA, 0.0015, 0.0120…
## $ cadmio                      <dbl> NA, NA, NA, NA, 0.003, 0.003, 0.001, 0.001…
## $ chumbo                      <dbl> NA, NA, NA, NA, NA, NA, NA, 0.002, 0.016, …
## $ cloreto                     <dbl> 27.80, 10.10, 9.60, 14.30, 10.10, 8.60, 8.…
## $ clorofila_a                 <dbl> NA, 2.19, 7.64, 40.10, 5.53, 3.56, 2.67, 1…
## $ cobre                       <dbl> NA, NA, NA, NA, 0.0040, 0.0020, 0.0030, 0.…
## $ coliformes_termotolerantes  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ coliformes_totais           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ condutividade               <dbl> 123.3, 81.0, 66.6, 139.1, 65.4, 69.0, 75.6…
## $ cromo_total                 <dbl> NA, NA, NA, NA, NA, NA, NA, 0.002, 0.009, …
## $ dbo                         <dbl> 3, 1, 2, 4, 2, 2, 2, 4, 4, 1, 3, 3, 1, 2, …
## $ demanda_quimica_de_oxigenio <dbl> 31, 15, 27, 24, 14, 17, 32, 40, 38, 23, 22…
## $ e_coli                      <dbl> 173.1, 28.8, 122.3, 770.1, 49.6, 34.7, 344…
## $ ferro                       <dbl> NA, NA, NA, NA, 3.120, 3.550, 1.940, 2.960…
## $ fitoplancton_cianobacterias <dbl> 171, 47, 1368, 129, 438, 68, 242, 2901, 18…
## $ fosfato_orto                <dbl> 0.071, 0.101, 0.089, 0.184, 0.109, 0.118, …
## $ fosforo_total               <dbl> 0.105, 0.202, 0.220, 0.233, 0.113, 0.184, …
## $ manganes                    <dbl> NA, NA, NA, NA, 0.051, 0.045, 0.066, 0.112…
## $ mercurio                    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ niquel                      <dbl> NA, NA, NA, NA, 0.0055, 0.0140, 0.0055, 0.…
## $ nitrato                     <dbl> NA, NA, 0.15, NA, NA, NA, NA, NA, NA, NA, …
## $ nitrogenio_amoniacal        <dbl> 0.032, 0.135, 0.530, 0.603, 0.174, 0.163, …
## $ nitrogenio_organico         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ nitrogenio_total_kjeldahl   <dbl> NA, 0.910, 1.350, 0.684, 1.120, NA, 1.907,…
## $ oxigenio_dissolvido         <dbl> 5.11, 6.39, 6.77, 3.25, 7.63, 7.50, 3.56, …
## $ ph                          <dbl> 6.48, 7.10, 6.21, 6.67, 6.73, 6.71, 6.27, …
## $ profundidade_coleta         <dbl> 0.88, 0.20, 0.20, 0.20, 0.20, 0.20, 0.20, …
## $ profundidade_total          <dbl> 2.40, 2.75, 2.75, 1.75, 2.15, 2.75, 3.00, …
## $ salinidade                  <dbl> NA, 0.04, NA, 0.07, 0.03, 0.03, 0.03, 0.07…
## $ solidos_dissolvidos_totais  <dbl> 100, 74, 123, 100, 87, 57, 86, 152, 20, 11…
## $ solidos_suspensos_totais    <dbl> 5.0, 15.0, 5.0, 13.0, 5.0, 21.0, 11.0, 15.…
## $ solidos_totais              <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ temperatura_da_agua         <dbl> 14.40, 17.43, 24.89, 23.30, 16.01, 20.36, …
## $ temperatura_do_ar           <dbl> 13.00, 16.00, 19.50, 23.00, 17.16, 21.30, …
## $ transparencia_da_agua       <dbl> 0.4, 0.2, 0.3, 0.3, 0.3, 0.2, 0.3, 0.3, 0.…
## $ turbidez                    <dbl> 23.00, 12.00, 49.10, 23.50, 15.50, 79.60, …
## $ vazao_recurso_hidrico       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ zinco                       <dbl> NA, NA, NA, NA, 0.0150, 0.0240, 0.0190, 0.…

# str(plan_litoral_medio) #outra maneira de visualizar a eSTRutura 

6 Requisitos pra gerar os gráficos

6.1 Definir o theme

Definindo uma cor de letra e de fundo padrão a ser adotada nos gráficos.

theme_grafs <- function(bg = "white", 
                        coloracao_letra = "black") {
  theme(
    plot.title = 
      element_text(
        hjust = 0.5,
        color = coloracao_letra,
        size = 19),
    
    axis.title.x = 
      # element_text(
      # color = coloracao_letra,
      # size = 15,
      # angle = 0,),
      element_blank(),
    axis.title.y = element_text(
      color = coloracao_letra,
      size = 15,
      angle = 90),
    
    axis.text.x = element_text(
      color = coloracao_letra,
      size = 17),
    axis.text.y = element_text(
      color = coloracao_letra,
      size = 17,
      angle = 0),
    
    strip.background = element_rect(fill = bg,
                                    linetype = 1,
                                    size = 0.5,
                                    color = "black"),
    strip.text = element_text(size = 17),
    panel.background = element_rect(fill = bg),
    plot.background = element_rect(fill = bg),
    plot.margin = margin(l = 5, r = 10,
                         b = 5, t = 5)
  )
}

6.2 Criar função para gerar boxplots com percentis 20 e 80

f <- function(x) {
  r <- quantile(x, probs = c(0.10, 0.20, 0.50, 0.80, 0.90))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  return(r)
}

7 Gráficos

7.1 Alterar nome e/ou ordem das estações

Atenção: Tomar muito cuidado com a ordem que estão dispostas as estações. Deve-se sempre respeitar o sentido nascente -> foz.

Conferir a ordem das estações através do mapa interativo, QGIS ou Google Earth Pro.

Elas acabam ficando fora de ordem por conta do código de ottobacias Para reordenar as estações, atente-se à função scale_x_discrete.

Neste exemplo iremos aprender a como mudar o nome das estações e reordená-las.

As estações estão na ordem certa, mas quero alterar o nome delas, como faço?

---

Utilizando o code chunk abaixo nós alteramos os nomes das estações para PM (ou Ponto de Monitoramento) seguido de um número. O PM1 representa o ponto mais próximo das nascentes, enquanto o PM6 o mais próximo da foz.

scale_x_discrete(
  limits = c(
    "87332500",
    "87420130",
    
    "87420150",
    "87420350",
    "87420500",
    "87510010"
  ),
  labels = c(
    "PM1", "PM2", "PM3", "PM4", "PM5", "PM6"
  )
)+

Perceba que iremos trocar a ordem das estações. Antes a estação 87332500 era a primeira a ser listada dentro da função limits = c(), ela era, portanto, o PM1. Agora o PM1 será a 87420130. Compare o código e o gráfico abaixo com o anterior.

scale_x_discrete(
  limits = c(
    "87420130",
    "87332500",
    
    "87420150",
    "87420350",
    "87420500",
    "87510010"
  ),
  labels = c(
    "PM1", "PM2", "PM3", "PM4", "PM5", "PM6"
  )
)+

“Ok mas… e se eu não quiser alterar o código ottobacias para PM1, PM2, só a ordem das estações, o que devo fazer?” Basta replicar a sequência correta das estações de limits = c() dentro da função labels = c().

   scale_x_discrete(
     limits = c(
       "87420130",
       "87332500",
       "87420150",
       "87420350",
       "87420500",
       "87510010"
     ),
     labels = c(
       "87420130",
       "87332500",
       "87420150",
       "87420350",
       "87420500",
       "87510010"
     )
   )+

Obs: Para a Bacia Hidrográfica do Rio Gravataí esse processo é necessário.

Sempre tomar cuidado com os limites do eixo y. Deixei como padrão o R buscar automaticamente qual o valor mínimo e máximo daquele parâmetro, mas o ideal é que se ajuste caso a caso.

Caso queira alterar o tamanho dos outliers, alterar o size nesse trecho dos códigos:

ggbeeswarm::geom_quasirandom(
     # grouponX = FALSE,
     size = 1.2, 
     alpha = .25,
     width = .07,
   )+

Todos os gráficos foram gerados utilizando como padrão os limites da Resolução nº 357/05 do CONAMA para ambientes lóticos. Caso queira alterar os limites para adequar aos ambientes lênticos, deve-se editar os ymin e ymax de cada retângulo (rect) do referido parâmetro.

Não esquecer que o R entende a casa decimal como ponto, e não como vírgula.

annotate("rect",
         xmin = -Inf, xmax = Inf,
         ymin = 13.3, ymax = Inf,
         alpha = 1,
         fill = "#ac5079")+ #>pior classe
annotate("rect",
         xmin = -Inf, xmax = Inf,
         ymin = 3.7, ymax = 13.3,
         alpha = 1,
         fill = "#fcf7ab")+ #classe 3
annotate("rect",
         xmin = -Inf, xmax = Inf,
         ymin = 0, ymax = 3.7,
         alpha = 1,
         fill = "#8dcdeb")+ #classe 1

7.2 Parâmetros de Qualidade da Água

7.2.1 Oxigênio Dissolvido

(graf_od <- plan_litoral_medio %>%
   ggplot(
     aes(
       x = cod_estacao,
       y = oxigenio_dissolvido,
     )
   )+
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = -Inf, ymax = 2,
            alpha = 1,
            fill = "#ac5079")+ #>pior classe
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 2, ymax = 4,
            alpha = 1,
            fill = "#eb5661")+ #classe 4
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 4, ymax = 5,
            alpha = 1,
            fill = "#fcf7ab")+ #classe 3
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 5, ymax = 6,
            alpha = 1,
            fill = "#70c18c")+ #classe 2
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin= 6, ymax = Inf,
            alpha = 1,
            fill = "#8dcdeb")+ #classe 1
   stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(
     title = "Oxigênio Dissolvido",
     x= NULL,
     y="mg/L"
   )+
   ggbeeswarm::geom_quasirandom(
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
   scale_y_continuous(
     expand = expansion(mult = c(0,0)),
     n.breaks = 11,
     # limits = c(-0.3,21)
     limits = c(
       min(plan_litoral_medio$oxigenio_dissolvido, na.rm = TRUE),
       max(plan_litoral_medio$oxigenio_dissolvido, na.rm = TRUE)+1)
   )+
   # scale_x_discrete(
   #   limits = c(
   #     "87332500",
   #     "87420130",
   #     "87420150",
   #     "87420350",
   #     "87420500",
   #     "87510010"
   #   ),
   #   labels = c(
   #     "PM1", "PM2", "PM3", "PM4", "PM5", "PM6"
   #   )
   # )+
   geom_smooth(
     method = "lm",
     se = FALSE, #se deixar TRUE gera o intervalo de confiança de 95%
     aes(group = 1),
     alpha = 0.5,
     na.rm = TRUE,
     size = 1
   )+
   theme_grafs()
)

Oxigênio Dissolvido

7.2.2 DBO

(graf_dbo <- ggplot(plan_litoral_medio,
                    aes(x = cod_estacao,
                   y = dbo))+
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 10, ymax = Inf,
            alpha = 1,
            fill = "#ac5079")+ #>pior classe
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 5, ymax = 10,
            alpha = 1,
            fill = "#fcf7ab")+ #classe 3
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 3, ymax = 5,
            alpha = 1,
            fill = "#70c18c")+ #classe 2
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 0, ymax = 3,
            alpha = 1,
            fill = "#8dcdeb")+ #classe 1
   stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(title = "Demanda Bioquímica de Oxigênio",
        x="Estação",
        y="mg/L")+
   ggbeeswarm::geom_quasirandom(
     # grouponX = FALSE,
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
   # scale_x_discrete(limits = c("87398500", 
   #                             "87398980", 
   #                             "87398900", 
   #                             "87398950", 
   #                             "87405500", 
   #                             "87406900", 
   #                             "87409900"),
   #                  labels = c("PM1", "PM2", "PM3", "PM4", "PM5", "PM6", "PM7")
   # )+
   scale_y_continuous(expand = expansion(mult = c(0.03,0.03)),
                      n.breaks = 8,
                      limits = c(
                        min(plan_litoral_medio$dbo, na.rm = TRUE),
                        10
                        # max(plan_litoral_medio$dbo, na.rm = TRUE)
                        ),
                      trans = "log10")+
   geom_smooth(method = "lm",
               se=FALSE, #se deixar TRUE gera o intervalo de confiança de 95%
               aes(group=1),
               alpha=.5,
               na.rm = TRUE,
               size = 1)+
   theme_grafs()
)

Demanda Bioquímica de Oxigênio

7.2.3 Escherichia coli

(graf_ecoli <- plan_litoral_medio %>% 
   ggplot(aes(cod_estacao,
              e_coli))+
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 3200, ymax = Inf,
            alpha = 1,
            fill = "#ac5079")+ #>pior classe
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 800, ymax = 3200,
            alpha = 1,
            fill = "#fcf7ab")+ #classe 3
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 160, ymax = 800,
            alpha = 1,
            fill = "#70c18c")+ #classe 2
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 0, ymax = 160,
            alpha = 1,
            fill = "#8dcdeb")+ #classe 1
   stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(title = "*Escherichia coli*",
        x="Estação",
        y="NMP/100mL")+
   scale_y_continuous(expand = expansion(mult = c(0.01, 0.01)),
                      # n.breaks = 9,
                      n.breaks = 6,
                      limits = c(min(plan_litoral_medio$e_coli, na.rm = TRUE),
                                 max(plan_litoral_medio$e_coli, na.rm = TRUE)),
                      trans = "log10",
                      labels = scales::number_format(accuracy = 1,
                                                     decimal.mark = ",",
                                                     big.mark = " "))+
   ggbeeswarm::geom_quasirandom(
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
 # scale_x_discrete(limits = c("87398500",
 #                             "87398980",
 #                             "87398900",
 #                             "87398950",
 #                             "87405500",
 #                             "87406900",
 #                             "87409900"),
 #                  labels = c("PM1", "PM2", "PM3", "PM4", "PM5", "PM6", "PM7")
 # )+
 geom_smooth(method = "lm",
 se=FALSE, #stardand error = desvio padrão -> se deixar TRUE gera o intervalo de confiança de 95%
             aes(group=1),
             alpha= 0.5, #transparencia de 50%
             na.rm = TRUE, #remover NAs
             size = 1)+
 theme_grafs()+
   theme(
     axis.text.y = element_text(
       angle = 90,
       # size=15,
       # face=2
     ),
     plot.title = 
       element_markdown(
         hjust = 0.5,
         color = "black",
         size = 19),
   )
)

Escherichia-coli

7.2.4 Fósforo Total

(graf_ptot <- ggplot(plan_litoral_medio,
                     aes(cod_estacao,
                         fosforo_total))+
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 0.15, ymax = Inf,
            alpha = 1,
            fill = "#ac5079")+ #>pior classe
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 0.1, ymax = 0.15,
            alpha = 1,
            fill = "#fcf7ab")+ #classe 3
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 0, ymax = 0.1,
            alpha = 1,
            fill = "#8dcdeb")+ #classe 1
   stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(title = "Fósforo total",
        x="Estação",
        y="mg/L")+
   scale_y_continuous(expand = expansion(mult = c(0.03,0.03)),
                      n.breaks = 8,
                      limits = c(min(plan_litoral_medio$fosforo_total, na.rm = TRUE),
                                 max(plan_litoral_medio$fosforo_total), na.rm = TRUE),
                      trans = "log10",
                      labels = scales::number_format(accuracy = .001,
                                                     decimal.mark = ",",
                                                     big.mark = " ")
   )+
   ggbeeswarm::geom_quasirandom(
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
   # scale_x_discrete(limits = c("87398500", 
   #                             "87398980", 
   #                             "87398900", 
   #                             "87398950", 
   #                             "87405500", 
   #                             "87406900", 
   #                             "87409900"),
   #                  labels = c("PM1", "PM2", "PM3", "PM4", "PM5", "PM6", "PM7")
   # )+
   geom_smooth(method = "lm",
               se=FALSE, #se deixar TRUE gera o intervalo de confiança de 95%
               aes(group=1),
               alpha=.5,
               na.rm = TRUE,
               size = 1)+
   theme_grafs()
)

Fósforo total

7.2.5 Nitrogênio amoniacal

(graf_namon <- ggplot(plan_litoral_medio,
                 aes(cod_estacao,
                     nitrogenio_amoniacal))+
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 13.3, ymax = Inf,
            alpha = 1,
            fill = "#ac5079")+ #>pior classe
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 3.7, ymax = 13.3,
            alpha = 1,
            fill = "#fcf7ab")+ #classe 3
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 0, ymax = 3.7,
            alpha = 1,
            fill = "#8dcdeb")+ #classe 1
   stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(title = "Nitrogênio amoniacal",
        x="Estação",
        y="mg/L")+
   scale_y_continuous(expand = expansion(mult = c(0.01, 0.05)),
                      n.breaks = 9,
                      limits = c(min(plan_litoral_medio$nitrogenio_amoniacal, na.rm = TRUE),
                                 15
                                 # max(plan_litoral_medio$nitrogenio_amoniacal, na.rm = TRUE)
                                 ),
                      trans = "log10",
                      labels = scales::number_format(accuracy = .001,
                                                     decimal.mark = ",",
                                                     big.mark = " "))+
   ggbeeswarm::geom_quasirandom(
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
   # scale_x_discrete(limits = c("87398500", 
   #                             "87398980", 
   #                             "87398900", 
   #                             "87398950", 
   #                             "87405500", 
   #                             "87406900", 
   #                             "87409900"),
   #                  labels = c("PM1", "PM2", "PM3", "PM4", "PM5", "PM6", "PM7")
   # )+
   geom_smooth(method = "lm",
               se=FALSE, #se deixar TRUE gera o intervalo de confiança de 95%
               aes(group=1),
               alpha=.5,
               na.rm = TRUE,
               size = 1)+
   theme_grafs()
)

Nitrogênio Amoniacal

7.2.6 Turbidez

(graf_turb <- ggplot(plan_litoral_medio,
                   aes(cod_estacao,
                       turbidez))+
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 100, ymax = Inf,
            alpha = 1,
            fill = "#ac5079")+ #>pior classe
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 40, ymax = 100,
            alpha = 1,
            fill = "#fcf7ab")+ #classe 3
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 0, ymax = 40,
            alpha = 1,
            fill = "#8dcdeb")+ #classe 1
   stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(title = "Turbidez",
        x="Estação",
        y="UNT")+
   scale_y_continuous(expand = expansion(mult = c(0.05, 0.05)),
                      n.breaks = 8,
                      limits = c(
                        # 1,
                        min(plan_litoral_medio$turbidez, na.rm = TRUE),
                        # 500
                        max(plan_litoral_medio$turbidez, na.rm = TRUE)
                      ),
                      trans = "log10",
                      labels = scales::number_format(accuracy = 1,
                                                     decimal.mark = ",",
                                                     big.mark = " "))+
    ggbeeswarm::geom_quasirandom(
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
   # scale_x_discrete(limits = c("87398500", 
   #                             "87398980", 
   #                             "87398900", 
   #                             "87398950", 
   #                             "87405500", 
   #                             "87406900", 
   #                             "87409900"),
   #                  labels = c("PM1", "PM2", "PM3", "PM4", "PM5", "PM6", "PM7")
   # )+
   geom_smooth(method = "lm",
               se=FALSE, #se deixar TRUE gera o intervalo de confiança de 95%
               aes(group=1),
               alpha=.5,
               na.rm = TRUE,
               size = 1)+
   theme_grafs()
)

turbidez

7.2.7 pH

(graf_pH <- ggplot(plan_litoral_medio,
                 aes(cod_estacao,
                     ph))+
   annotate("rect",
            xmin=-Inf,
            xmax=Inf,
            ymin=-Inf,
            ymax=6,
            alpha=1,
            fill="#eb5661")+ #classe 4
   annotate("rect",
            xmin=-Inf,
            xmax=Inf,
            ymin=9,
            ymax=Inf,
            alpha=1,
            fill="#eb5661")+ #classe 4
   annotate("rect",
            xmin=-Inf,
            xmax=Inf,
            ymin=6,
            ymax=9,
            alpha=1,
            fill="#8dcdeb")+ #classe 1
   stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(title = "pH",
        x="Estação",
        y="")+
   scale_y_continuous(expand = expansion(mult = c(0.01, 0.01)),
                      n.breaks = 8,
                      limits = c(4,11),
                      labels = scales::number_format(accuracy = 1,
                                                     decimal.mark = ",",
                                                     big.mark = " ")
                      )+
    ggbeeswarm::geom_quasirandom(
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
   # scale_x_discrete(limits = c("87398500", 
   #                             "87398980", 
   #                             "87398900", 
   #                             "87398950", 
   #                             "87405500", 
   #                             "87406900", 
   #                             "87409900"),
   #                  labels = c("PM1", "PM2", "PM3", "PM4", "PM5", "PM6", "PM7")
   # )+
   geom_smooth(method = "lm",
               se=FALSE, #se deixar TRUE gera o intervalo de confiança de 95%
               aes(group=1),
               alpha=.5,
               na.rm = TRUE,
               size = 1)+
   theme_grafs()
)

pH

7.2.8 Sólidos Totais

O parâmetro Sólidos Totais não teve medições no período analisado, então fica armazenado na tabela como valor vazio, ou rNA`. Quando você for fazer análise e tiver valores, altere o eval para TRUE.

(graf_solidos_totais <- ggplot(plan_litoral_medio,
                               aes(cod_estacao,
                                   solidos_totais))+
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 500, ymax = Inf,
            alpha = 1,
            fill="#ac5079")+ #>pior classe
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = -Inf, ymax = 500,
            alpha = 1,
            fill="#8dcdeb")+ #classe 1
stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(title = "Sólidos totais",
        x="Estação",
        y="")+
   scale_y_continuous(expand = expansion(mult = c(0.01, 0.05)),
                      n.breaks = 8,
                      limits = c(0,
                                 max(plan_litoral_medio$solidos_totais, na.rm = TRUE)
                                 ),
                      labels = scales::number_format(accuracy = 1,
                                                     decimal.mark = ",",
                                                     big.mark = " "))+
   ggbeeswarm::geom_quasirandom(
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
   # scale_x_discrete(limits = c("87398500",
   #                             "87398980",
   #                             "87398900",
   #                             "87398950",
   #                             "87405500",
   #                             "87406900",
   #                             "87409900"),
   #                  labels = c("PM1", "PM2", "PM3", "PM4", "PM5", "PM6", "PM7")
   # )+
   geom_smooth(method = "lm",
               se=FALSE, #se deixar TRUE gera o intervalo de confiança de 95%
               aes(group=1),
               alpha=.5,
               na.rm = TRUE,
               size = 1)+
   theme_grafs()
)

7.2.9 Condutividade

(graf_cond_elet <- ggplot(plan_litoral_medio,
                          aes(cod_estacao,
                              condutividade))+
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 500, ymax = Inf,
            alpha = 1,
            fill = "#eb5661")+ #classe 4
   annotate("rect",
            xmin = -Inf, xmax = Inf,
            ymin = 0, ymax = 500,
            alpha = 1,
            fill = "#8dcdeb")+ #classe 1
  stat_summary(
     fun.data = f,
     geom = 'errorbar',
     width = 0.3,
     position = position_dodge(width = 0.65),
   )+
   stat_summary(
     fun.data = f,
     geom = "boxplot",
     width = 0.7,
     fill = '#F8F8FF',
     color = "black",
     outlier.shape = 1, #se deixar NA fica só o jitter, se não, deixa 1
   )+
   # facet_wrap(~periodo)+
   labs(title = "Condutividade elétrica",
        x="Estação",
        y="µmhos/cm")+
   scale_y_continuous(expand = expansion(mult = c(0.05, 0.05)),
                      n.breaks = 8,
                      limits = c(min(plan_litoral_medio$condutividade, na.rm = TRUE),
                                 max(plan_litoral_medio$condutividade, na.rm = TRUE)),
                      trans = "log10",
                      labels = scales::number_format(accuracy = 1,
                                                     decimal.mark = ",",
                                                     big.mark = " "))+
   ggbeeswarm::geom_quasirandom(
     size = 1.2,
     alpha = .25,
     width = .07,
   )+
   # scale_x_discrete(limits = c("87398500", 
   #                             "87398980", 
   #                             "87398900", 
   #                             "87398950", 
   #                             "87405500", 
   #                             "87406900", 
   #                             "87409900"),
   #                  labels = c("PM1", "PM2", "PM3", "PM4", "PM5", "PM6", "PM7")
   # )+
   geom_smooth(method = "lm",
               se=FALSE, #se deixar TRUE gera o intervalo de confiança de 95%
               aes(group=1),
               alpha=.5,
               na.rm = TRUE,
               size = 1)+
   theme_grafs()+
   theme(
     axis.text.y = element_text(
       angle = 90,
       # size=15,
       # face=2
     )
   )
)

condutividade-eletrica

7.3 Análise ao longo do tempo

Para gerar um gráfico ao longo do tempo pra cada uma das estações é necessário alterar o ano_inicial e o ano_final. Também é necessário selecionar qual parâmetro que se quer fazer a visualização.

ano_inicial <- 2015
ano_final <- 2022

(timeline <- plan_litoral_medio %>%
  filter(ano_coleta > ano_inicial &
           ano_coleta <= ano_final) %>%
  dplyr::select(cod_estacao, e_coli, data_coleta) %>%
  group_by(cod_estacao) %>%
  ggplot(
    aes(x = data_coleta,
        y = e_coli,
        color = cod_estacao
    ))+
    geom_line(
      # aes(color = CODIGO),
      na.rm = TRUE)+
    geom_point(
      # aes(color = CODIGO),
      na.rm = TRUE)+
    scale_x_date(
      limits = as.Date(c(
        ymd(glue("{ano_inicial}-01-01")),
        ymd(glue("{ano_final}-01-01"))
        # NA #pode usar NA também
      )),
      expand = c(0.0, 0.0),
      date_breaks = "2 years",
      minor_breaks = "1 years",
      date_labels = "%Y",
    )+
  # geom_smooth(
  #   method = "lm", #regressao linear
  #   se = TRUE, #se deixar TRUE gera o intervalo de confiança de 95%
  #   aes(group = 1),
  #   alpha =.5,
  #   na.rm = TRUE,
  #   size = 0.3,
  #   # fullrange = TRUE,
  #   show.legend = TRUE
  # )+
  stat_smooth(
    geom = "smooth",
    # span = 0.2,
    se = TRUE, #se deixar TRUE gera o intervalo de confiança de 95%
    # aes(group = 1),
    # alpha =.5,
    na.rm = TRUE,
    # size = 0.3,
    fullrange = TRUE,
    show.legend = TRUE
  )+
  facet_wrap(
    ~cod_estacao,
    nrow = 4,
  )+
  theme_bw()
)

7.4 Correlação

parametros_iqa_lit_medio <- plan_litoral_medio %>%
  dplyr::select(cod_estacao,
         ph,
         dbo,
         e_coli,
         nitrogenio_amoniacal,
         # nitro_kjeldahl,
         # nitro_total,
         fosforo_total,
         temperatura_da_agua,
         turbidez,
         solidos_totais,
         oxigenio_dissolvido,
         condutividade,
         ano_coleta
         ) 

parametros_iqa_lit_medio %>% 
  dplyr::select(
    -cod_estacao,
    -ano_coleta,
    -solidos_totais
    ) %>%
  # group_by(cod_estacao) %>%
  rename(
    CE = condutividade,
    OD = oxigenio_dissolvido,
    # ST = solidos_totais,
    Turb = turbidez,
    Temp = temperatura_da_agua,
    Ptot = fosforo_total,
    NAmon = nitrogenio_amoniacal,
    pH = ph,
    DBO = dbo,
    E_coli = e_coli
    # NTK = nitro_kjeldahl
  ) %>% 
  ggcorr(
    method =
    "complete.obs",
    # "pearson",
    # "pairwise",
    name = "Correlação",
    label = TRUE,
    label_alpha = TRUE,
    digits = 3,
    low = "#3B9AB2",
    mid = "#EEEEEE",
    high = "#F21A00",
    # palette = "RdYlBu",
    layout.exp = 0,
    legend.position = "left",
    label_round = 3,
    # legend.size = 18,
    geom = "tile",
    nbreaks = 10,
  )+
  labs(title = "Correlação entre parâmetros físico-químicos na\nBacia Hidrográfica do Litoral Médio")+
  theme_linedraw()+
  theme(
    legend.position = c(0.15, 0.6),
    legend.title = element_text(size = 16),
    legend.text = element_text(size = 14),
    # legend.spacing = unit(element_text(),
                          # units = 5)
    plot.title = element_text(hjust = 0.5,
                              size = 16)
  )

correlação-parametros-qualidade-agua

7.4.1 Correlação entre parâmetros que compõem o IQA com significância estatística

Esse processo demanda bastante processamento, para desabilitá-lo deixe eval = FALSE.

correl_IQA_lit_medio <- parametros_iqa_lit_medio %>%
  dplyr::select(-cod_estacao) %>%
  ggpairs(title = "Correlação entre parâmetros que compõem o IQA",
          axisLabels = "show")
correl_IQA_lit_medio

plot

7.5 Salvar os gráficos

Replicar esse modelo pros gráficos que deseja salvar.
1. Criar novo code chunk (Ctrl+Alt+I)
2. Copiar a fórmula abaixo
3. Colar nesse novo code chunk

Os gráficos ficarão salvos em uma pasta que irá ser criada a partir do código.

ggsave("graf_od.png",
       plot = graf_od, #alteracao
       path = "./graficos",
       dpi = 300,
       type = "cairo")

8 Mapas

8.1 Pacotes pra geração de mapas

pacman::p_load(raster, leaflet, sf
               # prettymapr, rjson, rosm,
               # ggspatial
               # rgdal, rgeos,
               # gtools, tidyverse, rnaturalearth,
               # rnaturalearthdata, reticulate, maptools,
               # maps, ggplot2, ggspatial, rgeos, ggmap
               )

8.2 Mapa estático

Importando as informações necessárias pra gerar mapas de precipitação.
Fonte: GADM ¹, ²

Brasil <- getData(
  'GADM',
  country = 'Brazil',
  level = 3
) %>%
  st_as_sf()

## Warning in getData("GADM", country = "Brazil", level = 3): getData will be removed in a future version of raster
## . Please use the geodata package instead

RS <- subset(Brasil,
             NAME_1 == "Rio Grande do Sul")

lbl <- data.frame(month_abb = month.abb,
                  mes = 1:12)

plan_litoral_medio %>%
  ggplot(
    # aes(x = longtitude,
    #     y = latitude,
    #     map_id = region)
  )+
  geom_sf(
    data = RS
  )+
  theme_bw()

plot

Mapa estático do RS

# Definindo o SRC
RS <- RS %>%
  st_transform (crs = 4674) #4326 = WGS84, 4674 = SIRGAS2000

8.3 Mapa interativo com localização dos pontos de monitoramento

leaflet(RS) %>% 
  addProviderTiles(
    "Esri.WorldImagery" #Imagem de satélite
    # "OpenStreetMap.Mapnik" #OpenStreetMap -> Software livre
  ) %>% 
  addCircleMarkers(
    data = plan_litoral_medio,
    lng = ~longitude,
    lat = ~latitude,
    popup = ~paste(
      "<b>Estação:</b>",{cod_estacao},"<br>",
      "<b>Recurso hídrico:</b>", recurso_hidrico, "<br>",
      "<b>Município:</b>", municipio,
      sep = " "
    )
  ) 

Mapa interativo

8.4 Mapa de precipitação anual

Esse processo demanda bastante processamento, para desabilitá-lo deixe eval = FALSE.

st_centroid(RS)

output

## Simple feature collection with 1174 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -57.33034 ymin: -33.6509 xmax: -49.78182 ymax: -27.14727
## Geodetic CRS:  SIRGAS 2000
## First 10 features:
##       GID_0 NAME_0    GID_1            NAME_1 NL_NAME_1      GID_2     NAME_2
## 28115   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.1_1 Água Santa
## 28139   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.1_1 Água Santa
## 28163   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.1_1 Água Santa
## 28689   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.2_1      Agudo
## 28107   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.3_1  Ajuricaba
## 28119   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.3_1  Ajuricaba
## 27837   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.4_1    Alecrim
## 28607   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.5_1   Alegrete
## 28840   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.5_1   Alegrete
## 28842   BRA Brazil BRA.21_1 Rio Grande do Sul      <NA> BRA.21.5_1   Alegrete
##       NL_NAME_2        GID_3                         NAME_3 VARNAME_3 NL_NAME_3
## 28115      <NA> BRA.21.1.1_1                     Água Santa      <NA>      <NA>
## 28139      <NA> BRA.21.1.2_1                 Engenho Grande      <NA>      <NA>
## 28163      <NA> BRA.21.1.3_1 Santo Antonio dos Pinheirinhos      <NA>      <NA>
## 28689      <NA> BRA.21.2.1_1                          Agudo      <NA>      <NA>
## 28107      <NA> BRA.21.3.1_1                      Ajuricaba      <NA>      <NA>
## 28119      <NA> BRA.21.3.2_1                     Medianeira      <NA>      <NA>
## 27837      <NA> BRA.21.4.1_1                        Alecrim      <NA>      <NA>
## 28607      <NA> BRA.21.5.1_1                       Alegrete      <NA>      <NA>
## 28840      <NA> BRA.21.5.2_1                        Central      <NA>      <NA>
## 28842      <NA> BRA.21.5.3_1                          Leste      <NA>      <NA>
##         TYPE_3 ENGTYPE_3 CC_3 HASC_3                    geometry
## 28115 Distrito  District <NA>   <NA> POINT (-52.01191 -28.20432)
## 28139 Distrito  District <NA>   <NA>  POINT (-52.12523 -28.2061)
## 28163 Distrito  District <NA>   <NA> POINT (-52.08645 -28.26863)
## 28689 Distrito  District <NA>   <NA> POINT (-53.22639 -29.60797)
## 28107 Distrito  District <NA>   <NA> POINT (-53.77303 -28.21167)
## 28119 Distrito  District <NA>   <NA> POINT (-53.66801 -28.20777)
## 27837 Distrito  District <NA>   <NA> POINT (-54.78171 -27.65349)
## 28607 Distrito  District <NA>   <NA> POINT (-55.93205 -29.75749)
## 28840 Distrito  District <NA>   <NA> POINT (-55.78896 -29.78455)
## 28842 Distrito  District <NA>   <NA> POINT (-55.76203 -29.80277)

Prec <- getData(
  "worldclim",
  var = "prec",
  res = 0.5,
  lat = c(-30.033056, -29.68417), #procurar lat long no google e alterar aqui
  lon = c(-51.230000, -53.80694)
)

Prec_RS <- Prec %>%
  crop(RS) %>%
  mask(RS, na.rm = TRUE)

plot(Prec_RS)

plot

PPAnual_RS <- do.call("sum",
                       unstack(Prec_RS))
plot(PPAnual_RS)

plot

vls <- rasterToPoints(Prec_RS) %>%
  as_tibble() %>%
  gather(var, value, -x, -y) %>%
  mutate(mes = parse_number(var)) %>%
  inner_join(., lbl, by = 'mes') %>%
  dplyr::select(x, y, month_abb, value) %>%
  mutate(month_abb = factor(month_abb, levels = month.abb))

vls %>%
  filter(month_abb == 'Jan')

output

## # A tibble: 217,860 × 4
##        x     y month_abb value
##    <dbl> <dbl> <fct>     <dbl>
##  1 -53.1 -27.1 Jan         165
##  2 -53.0 -27.1 Jan         165
##  3 -53.0 -27.1 Jan         165
##  4 -53.0 -27.1 Jan         166
##  5 -53.0 -27.1 Jan         166
##  6 -53.4 -27.1 Jan         161
##  7 -53.4 -27.1 Jan         161
##  8 -53.4 -27.1 Jan         161
##  9 -53.1 -27.1 Jan         165
## 10 -53.1 -27.1 Jan         165
## # … with 217,850 more rows

summary(vls$value) #descobre o valor mínimo, médio e máximo de precipitação

output

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    39.0   128.0   141.0   140.3   154.0   206.0

8.5 Mapa com precipitação mensal

Caso queira gerar um mapa com a precipitação mensal, alternar o eval para TRUE. A geração desse mapa demanda bastante processamento, recomendável manter include = FALSE e echo = FALSE.

plot

9 Referências - Links úteis

Recomendo, abaixo, algumas referências. Também deixo uma trilha de aprendizado que julgo ser mais fácil do que a que percorri, conforme diagrama abaixo.

pacman::p_load(
  DiagrammeR
)
DiagrammeR::grViz(
  "digraph {
  graph [
  layout = dot,
  rankdir = LR
  ]

  node [shape = oval] #rectangle
  rec1 [label = 'Tentar reproduzir códigos\nde outras pessoas']
  rec2 [label = 'Entender a gramática\ndos gráficos']
  rec3 [label = 'Fazer um projeto\ndo zero']
  rec4 [label = 'Aprender RMarkdown']
  rec5 [label = 'procurar no Youtube/\nStackOverFlow']
  rec6 [label = 'pacotes']
  rec7 [label = 'tidyverse, ggplot2']
  rec8 [label = 'Desenvolver um Shiny App']

  # edge definitions with the node IDs
  rec1 -> rec5
  rec2 -> rec6 -> rec7
  rec1 -> rec2 -> rec3 -> rec4 -> rec8
  }",
  # height = 500
)

9.1 R

R for Data Science

Um dos principais livros para quem está querendo aprender R. Essa é a 1ª versão, publicada em 2017. Em 2023 será publicada a 2ª, através do link.

9.2 Cheatsheets

Cheatsheets

GitHub com as versões mais atualizadas

ggplot2 cheatsheet

ggplot2 cheatsheet em português

manipulação de dados em português

dplyr cheatsheet

9.3 Rmarkdown

Tutorial sobre RMarkdown

Rmarkdown basics

RMarkdown: The definitive guide

9.4 Pacotes

knitr

intro2r

• Getting started with R and RStudio

• Some R basics

• Data in R

• Graphics with ggplot

• Simple Statistics in R

• Reproducible reports with R markdown

9.5 Canais no Youtube

Pat Schloss #esse cara é genial

Ambiental Pro Dev

Marcelo Carvalho Ribeiro - Youtube

• Tem o Blog dele também.

Fernanda Peres

Marin Stats

Ainda precisa ser implementado/ajustado

• RSÁgua deve buscar dados anteriores a 2015.
• A Fepam tem os dados de altitude das estações. A partir disso pode ser calculado o % de saturação de Oxigênio Dissolvido, necessário para o cálculo do IQA.
• Tornar a coluna de chuva_24h como factor.
• A coluna hora_coleta está com formato incorreto.
• Gerar um Shiny Web App
• Tornar a coluna de municipio com apenas a primeira letra maiúscula.
• Sincronização via GitHub
• Aprender a gerar o mapa de precipitação mensal pra todo o RS
  • O RS está entre 2 fusos, então acaba “partindo no meio” a figura
• Futuramente: integração via Google Sheets.

Informações adicionais

E-mail para contato: leonardofwink@gmail.com

GitHub: https://github.com/leonardofwink/

LinkedIn: https://www.linkedin.com/in/leonardofwink/

---

1. https://gadm.org/maps/BRA.html↩︎

2. https://gadm.org/data.html↩︎