Análise Hidrológica com R: uma introdução para iniciantes: Manipulação dos dados

Pacotes

Nessa seção vamos usar mecanismos de pacotes já feitos no R. Existem inúmeros pacotes, cada um com suas funções e/ou bancos de dados para diferentes assuntos. Eles podem, inclusive, ser apenas um aglomerado de diferentes pacotes para não precisar baixar todos um por um. A primeira etapa de utilizar um pacote é a sua instalação, isso é feito com o comando install.packages() com o nome do pacote estando entre aspas (““) dentro dos parêntesis. Assim podemos fazer com todos que se encontram no CRAN (repositório que contém o que queremos). Na aula de hoje utilizaremos apenas o tidyverse, que é um agregado de pacotes extremamente populares hoje para manipulação e análise de dados.

install.packages("tidyverse")

Para que o R não fique muito pesado, todos os pacotes baixados precisam ser carregados na sua seção do R - isso é feito com o comando library() com o pacote estando detro do parêntesis sem aspas.

Podemos notar que ao carregar o pacote dentro da nossa seção primeiro vemos quais outros pacotes que foram carregados e suas respectivas versões. Nesse exemplo, ao carregar o tidyverse, ele contém o ggplot2, purrr, tibble, dplyr, tidyr, stringr, readr e forcats. Além disso, aparece uma mensagem mostrando “conflitos”. Esses conflitos são quando diferentes pacotes utilizam funções de mesmo nome (as vezes para fazer coisas diferentes), assim o R te avisa quais funções são iguais e quais serão as utilizadas normalmente no programa. É possível usar qualquer função se deixarmos explícito de qual pacote queremos (por exemplo stats::filter())

library(tidyverse)

Puxar arquivo anterior e analisar ele

Primeiro vamos baixar apenas uma das estações que fizemos na última aula. Ao final podemos repetir o que foi feito anteriormente e criar um loop para repetir para as outras 3 estações.

dados_ANA <- read.table(file = "dados/60435000.txt",
                        sep = "\t",
                        dec = ".",
                        header = T,
                        fileEncoding = "UTF-8")

Apesar de termos mexido nesses dados na última aula, o R não sabe disso, ele baixa o arquivo como se fosse um outro qualquer, assim a coluna da estação é do tipo “inteiro” e a segunda coluna (de data) está do tipo “caractere”. Vemos isso abaixo e em seguida corrigimos ambas para caractere e data, respectivamente.

class(dados_ANA$Cod_estacao)

[1] "integer"

class(dados_ANA$Data)

[1] "character"

dados_ANA$Cod_estacao <- as.character(dados_ANA$Cod_estacao)
dados_ANA$Data <- as.Date(dados_ANA$Data)

Manipulações básicas nos dados

Agora que temos um arquivo do jeito que queríamos (diferente do que pegamos diretamente do servidor da ANA), fazer algumas análises podem ser bem fáceis. Se quisermos saber de quando a quando essa estação tem dados, basta vermos os valores mínimos e máximos da Data.

min(dados_ANA$Data)

[1] "1978-05-01"

max(dados_ANA$Data)

[1] "2021-11-30"

Se quisermos saber qual valor mínimo, máximo e médio de vazão podemos também, mas aparece um problema! Nesses valores temos dados em branco, ou NA (Not Available). Dados faltantes sempre são problemático e devem ser tratados quando possível. Para algumas funções nativas, como max(), min() e mean() temos a opção de apenas ignorar eles com o argumento na.rm = TRUE, assim os valores máximos, médios e mínimos são calculados sem considerar os dados em branco.

max(dados_ANA$Vazao, na.rm = T)

[1] 37.4633

min(dados_ANA$Vazao, na.rm = T)

[1] 0.0539

mean(dados_ANA$Vazao, na.rm = T)

[1] 2.096619

Dados em Branco (NA)

Antes mesmo de começar a fazer essas avaliações é interessante e importante entender melhor os dados em branco. Primeiro, quantos existem? Para isso podemos usar a função is.na(), vamos rodar ela apenas para os primeiros 15 dados de vazão para mostrar na tela.

is.na(dados_ANA$Vazao)[1:15]

 [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[12] FALSE FALSE FALSE FALSE

Essa função retorna um dado TRUE ou FALSE para cada valor avaliado. Como dados lógicos também são binários/numéricos (TRUE = 1 e FALSE = 0), podemos somar esses valores para ver quantos dados em branco temos no total. Temos então sum(is.na(dados_ANA[,3])) dados em branco.

sum(is.na(dados_ANA$Vazao))

[1] 87

Para saber quais as datas que tem o valor de vazão em branco basta aplicarmos esse is.na() na coluna de datas.

dados_ANA$Data[is.na(dados_ANA$Vazao)]

 [1] "1978-05-01" "1978-05-02" "1978-05-03" "1978-05-04" "1978-05-05"
 [6] "1978-05-06" "1978-05-07" "1978-05-08" "1978-05-09" "1978-05-10"
[11] "1978-05-11" "2008-12-01" "2008-12-02" "2008-12-03" "2008-12-04"
[16] "2008-12-05" "2008-12-06" "2008-12-07" "2008-12-08" "2008-12-09"
[21] "2008-12-10" "2008-12-11" "2008-12-12" "2008-12-13" "2008-12-14"
[26] "2008-12-15" "2008-12-16" "2008-12-17" "2008-12-18" "2008-12-19"
[31] "2008-12-20" "2008-12-21" "2008-12-22" "2008-12-23" "2008-12-25"
[36] "2008-12-26" "2008-12-27" "2008-12-28" "2008-12-29" "2008-12-30"
[41] "2008-12-31" "2009-01-01" "2009-01-02" "2009-01-03" "2009-01-04"
[46] "2009-01-05" "2009-01-06" "2009-01-07" "2009-01-08" "2009-01-09"
[51] "2009-01-10" "2009-01-11" "2009-01-12" "2009-01-13" "2009-01-14"
[56] "2009-01-15" "2019-08-01" "2019-08-02" "2019-08-03" "2019-08-04"
[61] "2019-08-05" "2019-08-06" "2019-08-07" "2019-08-08" "2019-08-09"
[66] "2019-08-10" "2019-08-11" "2019-08-12" "2019-08-13" "2019-08-14"
[71] "2019-08-15" "2019-08-16" "2019-08-17" "2019-08-18" "2019-08-19"
[76] "2019-08-20" "2019-08-21" "2019-08-22" "2019-08-23" "2019-08-24"
[81] "2019-08-25" "2019-08-26" "2019-08-27" "2019-08-28" "2019-08-29"
[86] "2019-08-30" "2019-08-31"

E se quiser o último (a última data com vazão em branco), podemos fazer o max() dessas datas todas.

max(dados_ANA$Data[is.na(dados_ANA$Vazao)])

[1] "2019-08-31"

O que fazer com cada dado em branco vai depender do que se quer com os dados. Como completar essas falhas é uma temática extremamente relevante e difícil de responder até mesmo hoje. Pararemos por hora de falar desses valores, e não entraremos muito a fundo sobre o que fazer com eles (até porque depende do que se quer fazer) nesse curso.

Ano hidrológico

Em hidrologia estatística, é muito comum utilizarmos o ano hidrológico ao invés do calendário “normal” (janeiro a dezembro) para fazer análises estatísticas. Onde começa o ano hidrológico também não é tão simples, mas a via de regra é: entre as secas e cheias.

Faremos aqui uma função para transformar nossa data em ano hidrológico, assim conseguiremos fazer algumas análises anuais interessantes. A função, chamada aqui de fun_ano_hidro() tomará 2 argumentos: um sendo as datas que queremos transformar e o outro em que mês que se inicia o ano hidrológico, aqui deixamos pré-definido como 8, mas se quiser mudar isso, basta alterar o argumento na hora de chamar ela (assim não ficamos preso a esse número)! Primeiro devemos garantir que as datas estão no formato certo, portanto fazemos essa transformação no começo da função. Em seguida definimos qual será o ano hidrológico. O que fazemos é avaliar se o mês em questão é menor ou maior do que o começo do ano hidrológico definido pela função. Se for maior, significa que já entramos em outro ano hidrológico, portanto o ano em questão é o seguinte (o mês 8, Agosto, de 2000 entra para o ano hidrológico de 2001 e assim por diante).

fun_ano_hidro <- function(datas,
                          comeco_ano_hidro = 8){
 
  # Se não tiver no formato de datas, fazer essa correção
  datas <- as.Date(datas)
  
  # Se vai pular (deslocar) um ano ou não
  desloc_ano <- ifelse(month(datas) < comeco_ano_hidro, 0, 1)
  
  # Ano Hidro
  ano_hidro <- year(datas) + desloc_ano
  
  return(ano_hidro)
}

Abaixo testamos ela fazendo uma nova coluna no nosso dataframe.

dados_ANA$ano_hidro <- fun_ano_hidro(dados_ANA$Data)

head(dados_ANA)

  Cod_estacao       Data Vazao ano_hidro
1    60435000 1978-05-01    NA      1978
2    60435000 1978-05-02    NA      1978
3    60435000 1978-05-03    NA      1978
4    60435000 1978-05-04    NA      1978
5    60435000 1978-05-05    NA      1978
6    60435000 1978-05-06    NA      1978

Tabela Resumo

Aqui introduzimos um novo conceito do pacote magrittr (ele está internamente no tiyverse). Esse conceito se chama pipe operator e é dado pelo símbolo %>%. Basicamente o que ele faz é tomar o que está na esquerda dele, e usar como argumento da função a direita. Assim ao invés de escrever as funções umas dentro das outras, podemos deixar o código mais limpo e fácil de entender. Olhe o exemplo abaixo. Nele estamos tirando o log das vazões, com a função log(), depois calculando as diferenças desse log em um passo de tempo com diff(), calculamos o exponencial disso com exp(), em seguimas calculamos a média desses valores e depois arredondamos isso para ter apenas 2 casas decimais.

round(mean(exp(diff(log(dados_ANA$Vazao))), na.rm = T), 2)

[1] 1.03

É fácil se perder lendo isso pois no R normalmente lemos de dentro para fora das funções (e portanto da direita para a esquerda) e é mais fácil ainda confundir qual argumento entra em qual função. Usando o pipe operator, dado por %>%, o código fica mais facilmente interpretável. Selecionamos a vazão e em seguida tiramos o log, depois a diferença, depois o expoente… Quando uma função tem diversos argumentos podemos escrever explicitamente onde que o objeto a esquerda do pipe vai entrar colocando um ponto no lugar dele, não precisando que ele sempre entre no primeiro argumento da função.

dados_ANA$Vazao %>%
  log() %>%
  diff() %>%
  exp() %>%
  mean(., na.rm = T) %>%
  round(digits = 2, x = .)

[1] 1.03

Agora que temos o ano hidrológico para cada data, podemos fazer a nossa série de máximos, mínimos e médias anuais a partir dos dados diários. Primeiro cria-se um dataframe vazio com as colunas para cada ano, uma coluna para os máximos, outra para os mínimos e outra para as médias anuais. Podemos também criar uma quarta coluna para o número de valores NA nesse período.

Já iremos fazer isso no formato de função, para que possamos replicar para as outras estações também! A função tomará como argumento a tabela original de dados (dados_ANA), e irá rodar a função fun_ano_hidro() internamente. Assim teremos uma nova coluna criada com a função mutate(). Em seguida agrupamos os dados por estação (se tivesse mais de uma por exemplo) e por ano hidrológico usando o group_by e falamos o que faremos ao juntar esses dados com o summarise(). Um problema é se o ano hidrológico estiver todo com falha, nesses casos as funções usadas retornarão Inf (infinito), -Inf ou NaN (Not a Number), portanto depois de agruparmos, iremos trocar qualquer valor não-finito na tabela por NA.

Rodar uma outra função internamente requer alguns cuidados. Por exemplo, aqui assumimos que a tabela virá com uma coluna chamada “Data”. Caso isso não aconteca a nossa função não irá funcionar.

fun_resumo_ano <- function(dados_hidro){
  
  tabela_resumo <-
    dados_hidro %>%
    mutate(ano_hidro = fun_ano_hidro(Data)) %>%
    group_by(Cod_estacao, ano_hidro) %>%
    summarise(maxima = max(Vazao, na.rm = T),
              minima = min(Vazao, na.rm = T),
              media = mean(Vazao, na.rm = T),
              NAs = sum(is.na(Vazao)))
  
  tabela_resumo[!is.finite(tabela_resumo$maxima), c(3, 4, 5)] <- NA
  
  return(tabela_resumo)
}

Tabela final (loop)

Agora podemos aplicar essa função diretamente aos dados_ANA da última aula. Na realidade, podemos fazer um loop para abrir cada uma das 4 estações, rodar essa função e juntar essas tabelas em uma só com a função rbind() - essa função junta diferentes dataframes pelas mesmas colunas, colocando um embaixo do outro. Abaixo rodamos o loop de outra maneira, definindo antes quais os códigos das estações e rodando o for() variando de 1 a 4. Primeiro lemos a primeira estação e, se for a primeira iteração, salvamos ela em um objeto chamado tabela_final. Se não for a primeira iteração, juntamos elas com o rbind().

estacoes_cod <- c(60435000, 60436000, 60436190, 60443000)

for(i in 1:4){
  dados_proxy <- read.table(file = paste0("dados/", estacoes_cod[i],".txt"),
                            sep = "\t",
                            dec = ".",
                            header = T,
                            fileEncoding = "UTF-8")
  
  dados_proxy <- fun_resumo_ano(dados_proxy)
  
  if(i == 1) tabela_final <- dados_proxy
  if(i != 1) tabela_final <- rbind(tabela_final, dados_proxy)
}

Por último, podemos salvar esse arquivo também, para usarmos nas aulas seguintes.

write.table(x = tabela_final,
            file = "Resumo_estacoes.txt",
            sep = "\t",
            dec = ".",
            row.names = FALSE,
            fileEncoding = "UTF-8")

Manipulação dos dados