OBJETIVOS

Os principais objetivos deste post são:

CARREGANDO OS PACOTES

Antes de começarmos a brincadeira, vamos importar e carregar os pacotes que serão necessários. Faço uma breve descrição deles a seguir, quem tiver com pressa, é só carregar carregar os pacotes e seguir em frente.

Bem, agora que tivemos uma ideia (ainda que meio vaga, mas suficiente para o que precisamos) dos pacotes que vamos usar, vamos carregá-los. Abaixo o código em R para carregar os pacotes. Quem ainda não tem esses pacotes instalados, é preciso instalá-los primeiro antes de carregá-los.

##########################
##Instalando os pacotes## 
##########################
#Para quem usa o Rstudio, pode apenas ir na "aba" packages, clilar em install e digitar o nome dos pacotes um a um, se preferir). Eu deixei em comentário as linhas abaixo, quem quiser instalar diretamente pelo script baixa tirar o comentário (#) da frente do texto.

#install.packages("ggplot2",repos = "http://cran.us.r-project.org")

#install.packages("dplyr",repos = "http://cran.us.r-project.org")

#install.packages("lubridate",repos = "http://cran.us.r-project.org")

#install.packages("reshape2",repos = "http://cran.us.r-project.org")

#install.packages("gganimate",repos = "http://cran.us.r-project.org")

##########################
##Carregando os pacotes## 
##########################

#Depois de instalados precisamos carregar os pacotes, caso eles não estejam carregados ainda

library(dplyr)
library(lubridate)
library(reshape2)
library(ggplot2)
library(gganimate)

LENDO E MANIPULANDO OS DADOS DE ENTRADA

Baixando os Dados

Antes de iniciarmos a manipulação dos dados, é preciso acessar o arquivo com os dados. Os dados consistem em um histórico diário de Evapotranspiração Potencial (E), de Precipitação (P) e de vazão observada (Q). Eu deixei os dados para download em uma pasta do meu repositório na plataforma GitHub. Os dados podem ser baixados clicando neste link:https://github.com/cassiorampinelli/postsData/blob/master/EPQd.txt . E na sequência clicando no botão “Raw”. Uma nova página html irá abrir com o conteúdo do arquivo. Basta clicar na página com o botão direito e salvar no diretório onde o seu arquivo de trabalho do R está localizado. Salve com o nome “EPQd.txt”.

Definindo a pasta de trabalho

Salve o seu script do R em uma pasta e defina esta pasta como o seu diretório corrente de trabalho. Para isso basta usar o comando setwd e inserir o endereço da pasta nele. Certifique-se que consiste no mesmo diretório onde você salvou o arquivo EPQd.txt. No meu caso eu vou “setar” como local de trabalho uma pasta local. Basta você substituir o seu endereço local no código abaixo.(Observe o padrão das barras que normalmente é ao contrário do que o windows utiliza. Acho que no Mac é oposto)

setwd("E:/Users/Cassio/Docs/posts/2019_03_11")

Carregando e olhando os dados

Agora que salvamos os dados na nossa pasta de trabalho vamos carregá-lo com o comando read.table e dar uma olhada nos dados com o comando head que mostra algumas linhas iniciais da tabela. Veja que a primeira coluna contem os dias em que houve observação. A segunda coluna contém a série de evapotranspiração potencial (etp), a terceira coluna a série de chuva (pluv) e a última coluna as vazões (discharge). As unidades da evapotranspiração potencial e chuva são mm/dia e das vazões m³/s.

#Lendo os dados e mantendo os cabeçalhos
EPQOd<-read.table(file="EPQd.txt",header=T)

#Dando uma olhada no que tem nas primeiras linhas da tabela
head(EPQOd)
##         date      etp         plu discharge
## 1 1983-05-01 3.224745 0.009977963   96.5509
## 2 1983-05-02 3.179469 0.019636667   94.9420
## 3 1983-05-03 2.778005 0.013762269   94.1411
## 4 1983-05-04 1.785748 0.000000000   91.7529
## 5 1983-05-05 2.347023 0.084731326   90.9616
## 6 1983-05-06 2.879173 0.003784305   89.3864

Convertendo a série diária em série mensal

Bem, se olharmos quantas linhas tem esta tabela veremos que temos um histórico muito grande de dados diários iniciando em 01/05/1983 até 31/12/2012 (com algumas falhas, valores com NA indicam falhas).

Mas os dados estão em passos de tempo diário, como convertemos para passos de tempo mensais, sem termos que fazer isto na “unha”?

É aí que entram as funções dos pacotes dplyr e lubridate. Vamos fazer o seguinte:

  • Definir nosso intervalo de dias por meio de uma variável que vamos chamar period;

  • Para o histórico diário de evapotranspiração e chuva, a cada mês vamos somar os valores para ter o valor mensal acumulado. Por exemplo, a precipitação do mês de maio/1983 será a soma de todos os valores diários de precipitação para o mês de maio, do mês de setembro/1983, a soma de todos os valores diários de precipitação do mês de setembro, e assim sucessivamente. Para as vazões, vamos calcular a média das vazões diárias para cada mês e agrupá-los.

  • Converter o formato de dados de data frame para matriz

Abaixo o código para fazer essa “mágica” com alguns comentários nas linhas de código.

#Definindo todo o período de dados
period<-data_frame(date=seq(as.Date("1983-05-01"), as.Date("2012-12-31"), by=1))

#Trocando a primeira coluna do arquivo EPQOd pela série diária no formato period
EPQ=cbind(period,EPQOd[,2],EPQOd[,3],EPQOd[,4])

#Renomeando as colunas
colnames(EPQ) <- c("date", "etp","plu","discharge")

#Reagrupando os dados diários em mensais
EPQm=EPQ %>% group_by(month=floor_date(date, "month"))%>%summarize(etp=sum(etp),plu=sum(plu),discharge=mean(discharge))

#Salvando os dados organizados mensalmente
EPQOm=EPQm

Agora vamos olhar como ficou a série reagrupada na variável EPQOm, veja que agora cada linha apresenta os valores acumulados mensais para etp e plu e as médias para discharge. (Obs, a sintaxe usa o formato 1983-05-01 para indicar o mês de maio/1983, 1983-06-01 para o mês de junho, e assim sucessivamente).

Uma mão na roda esta função eim?!! Ao menos melhor do que fazer isso na “unha” ou usando o excel…

#Olhando as primeiars linhas dos dados agrupados mensalmente
head(EPQOm)
## # A tibble: 6 x 4
##   month        etp    plu discharge
##   <date>     <dbl>  <dbl>     <dbl>
## 1 1983-05-01  72.6  63.4       88.7
## 2 1983-06-01  68.1  11.3       78.0
## 3 1983-07-01  73.9  35.0       66.5
## 4 1983-08-01  97.3   2.54      55.6
## 5 1983-09-01  86.2 103.        61.9
## 6 1983-10-01 105.  174.        91.6

Bem, agora vamos deixar os dados preparados para usarmos mais a frente. Vamos trabalhar apenas com os meses de maio/1983 até Dezembro/1987. Ou seja, 56 meses de dados. Criaremos uma matriz EPQ, contendo apenas os primeiros 56 meses de dados que nos interessam e removeremos a primeira coluna que indica os períodos.

#Removendo a primeira coluna
EPQ<-EPQOm[,-1]

#Selecionando o período de interesse
EPQ<-EPQ[1:56,]

#Convertendo a estrutura de data frame para matriz
EPQ<-as.matrix(EPQ)

#Checando os dados
head(EPQ)
##            etp        plu discharge
## [1,]  72.61226  63.374100  88.70539
## [2,]  68.12304  11.261616  78.00904
## [3,]  73.89580  35.003316  66.49162
## [4,]  97.31646   2.542564  55.62024
## [5,]  86.19283 102.528167  61.90446
## [6,] 105.00501 174.157001  91.61406

O MODELO HIDROLÓGICO SMAP MENSAL

O modelo hidrológico SMAP (Soil Moisture Accounting Procedure) proposto por Lopes et al. (1981) para passos de tempo mensais é um modelo do tipo chuva vazão conceitual que procura representar o armazenamento e os fluxos de água na bacia através de reservatórios lineares fictícios. O modelo possui quatro parâmetros (SAT, PES, CREC, K ) que usualmente são calibrados e dois parâmetros de iniciliazação (Tuin e Ebin).

A estrutura do modelo é composta por dois reservatórios que buscam representar o armazenamento e os fluxos na camada superior do solo e no aquífero, como pode ser visualizado na Figura 1, a seguir.

Figure 1: Modelo SMAP para passos de tempo mensais

Figure 1: Modelo SMAP para passos de tempo mensais

A cada evento de precipitação (P), realiza-se um balanço de massa na bacia hidrográfica em estudo. Uma parcela de (P) é transferida como escoamento superficial (ES), estimado por uma equação exponencial que depende de P, da taxa de umidade do solo (TU) e de um expoente (PES).

A lâmina restante da precipitação, subtraída do escoamento superficial (P-ES), sofre perda por evaporação (P-ES-EP) e é então adicionada a um reservatório, o qual representa a camada superior do solo. Neste reservatório, a umidade do solo é atualizada ao longo do tempo através das perdas por evapotranspiração real (ER), que dependem do nível do reservatório (RSOLO) e da capacidade de saturação do solo (SAT).

A saída deste segundo reservatório representa a recarga do reservatório subterrâneo (REC), estimada com base nataxa de umidade (TU), no nível do reservatório(RSOLO) e no coeficiente de recarga (CREC).

O nível d’água no segundo reservatório(RSUB)é então deplecionado a uma taxa constante de recessão do escoamento de base (K), resultando no escoamento de base (EB) propriamente dito. A soma de ES com EB fornece a vazão total no ponto de controle da bacia.

O modelo possui quatro parâmetros, a saber: capacidade de saturação do solo (SAT); expoente relacionado à geração de escoamento superficial (PES); o coeficiente de recarga do aquífero (CREC), que está relacionado com a permeabilidade da zona não saturada do solo; a taxa de deplecionamento (K) do nível d’água do segundo reservatório (RSUB), responsável pela geração do escoamento de base (EB).

Além disso, duas condições iniciais precisam ser definidas, a taxa de umidade inicial do solo (TUin), que determina o nível inicial do segundo reservatório (RSOLO) e o escoamento de base inicial (EBin).

Como dados de entrada são necessários um histórico de precipitação média mensal (em mm de chuva), evapotranspiração potencial mensal (em mm) e a área de drenagem da bacia hidrográfica modelada (em km²).

As linhas de código, a seguir, apresentam o passo a passo do modelo. Neste exemplo, criaremos uma função chamada SmapM, cujos dados de entrada serão uma matriz EPQ, contendo o histórico de evapotranspiração potencial, precipitação e vazão, um vetor para contendo os valores para os 4 parâmetros (SAT, PES, CREC, K) e Area uma variável numérica que contém a área de drenagem da bacia.

A saída do modelo será um objeto do tipo list, uma lista na sintaxe do R que contem 6 vetores com o histórico de vazões simuladas (Qcalc), os níveis de água no reservatório superficial Rsol(Rsol), no reservatório subterrâneo Rsub(rsub), os níveis equivalentes à evapotranspiração real (er) e as vazões correspondentes ao escoamento superficial (es) e escoamento de base (eb).

Abaixo as linhas de código comentadas (Ficou um pouco grande, mas a ideia é tentar comentar cada linha de código para ficar mais fácil de compreender o que significam para quem está tendo contato com o modelo pela primeira vez)

smapM<-function(EPQ,para,Area) 
  
{
  ###################################
  #Inicializando os dados de entrada 
  ####################################

  E<-EPQ[,1]#Cria um vetor "E" com o histórico de evapotranspiração potencial
  Prec<-EPQ[,2]#Cria um vetor "E" com o histórico de evapotranspiração potencial 
  n<-nrow(EPQ) #Número total de passos de tempo
  
  ###################################
  #Inicializando os parâmetros do modelo 
  ####################################
  
  #sat, pes, crec, k,tuin, ebin
  sat=para[1]#Armazena o elemento 1 do vetor para na variável sat 
  pes=para[2]#Armazena o elemento 2 do vetor para na variável pes
  crec=para[3] #Armazena o elemento 3 do vetor para na variável crec
  crecp=crec/100 #Transforma crec de percentual para centesimal
  k=para[4]  #Armazena o elemento 4 do vetor para na variável k   
  ke=0.5^(1/k)#Calcula a constante de recessão do escoamento superficia
  
  #Para este exemplo, vamos assumir os seguintes valores iniciais para umidade e escoamento de base:
  tuinp=0.4
  ebin=20
  
  ###################################
  #Inicializando as variáveis de estado 
  ####################################
  #OBS: Estas variáveis mudam a cada passo de tempo
  
  Rsol<-0 #Nível do reservatório Rsol
  rsub<-0 #Nível do reservatório Rsub
  tu<-0   #Taxa de umidade 
  dsol<-0 #Variação de nível do reservatório
  Qcalc<-0 #Vazão simulada
  es<-0    #Escoamento superficial
  er<-0   # Evapotranspiração real
  eb<-0   #Escoamento de base
  rec<-0  # Recarga do reservatório Rsub
  
  #Nível inicial dos reservatórios Rsol e Rsub
  Rsol[1]=tuinp*sat
  rsub[1]=ebin/(1-ke)/Area*2630
  
  #Taxa de umidade inicial
  tu[1]=tuinp
  
   i<-2 # Variável i é inicializada para o segundo passo de tempo
   
  ###################################
  #Loop para rodar o modelo para cada passo de tempo 
  ####################################

   while (i<=n+1) #Loop while que é processado até que o contador i atinga n+1

  {
    # Calcula a variação de nível do reservatorio do solo
    dsol[i]=0.5*(Prec[i-1]-Prec[i-1]*(tu[i-1]^pes)-E[i-1]*tu[i-1]-Rsol[i-1]*crecp*(tu[i-1]^4))
    
    #Calcula a taxa de umidade no passo i
    tu[i]=(Rsol[i-1]+dsol[i])/sat
    
    #Calcula o escoamento de base e verifica se há valores     não numéricos, fazendo as devidas correções.
    if(is.na(tu[i])){tu[i]=0}
    if(tu[i]<0){tu[i]=0}
    if(tu[i]>1){tu[i]=1}
    es[i]=Prec[i-1]*tu[i]^pes
    
    #Calcula a evaporação real
    er[i]=E[i-1]*tu[i]
    
    #Calcula a recarga do reservatorio subterrâneo
    rec[i]=Rsol[i-1]*crecp*(tu[i]^4)
    
    #Calcula o escoamento de base
    eb[i]=rsub[i-1]*(1-ke)
    
    #Atualiza o nível do     reservatório subterrâneo
    rsub[i]=rsub[i-1]-eb[i]+rec[i]
    
    #Atualiza o nível do reservatorio do solo e verifica      se o     mesmo é superior a capacidade do servatório      (sat), se há      valor não numérico (NA) ou negativo,     fazendo os ajustes       necessários.
    Rsol[i]=Rsol[i-1]+Prec[i-1]-es[i]-er[i]-rec[i]

      if(is.na(Rsol[i])){Rsol[i]=0}
      if(Rsol[i]>sat)
        {es[i]=es[i]+Rsol[i]-sat
        Rsol[i]=sat}
      if(Rsol[i]<0){Rsol[i]=0}

    #Calcula vazão total simulada
    Qcalc[i-1]=(es[i]+eb[i])*Area/2630

     i<-i+1 #atualiza o contador
    }  

  #Cria uma lista com os valores calculados
  output=list(Qcalc,Rsol,rsub,er,es,eb)
  
  return(output) # Retorna a lista com todos os valores calculados
  
}

RODANDO O SMAP

Bem, agora que já temos a função do nosso modelo SMAP e, também, os dados de entrada. Podemos rodar o modelo, chamando a função que configuramos no passo anterior e passando os argumentos requeridos pela função.

Vamos assumir que nossa área de drenagem seja 4858 km² e vamos adotar como valores para os parâmetros do modelo os seguintes os seguintes:

  sat=1795.2

  pes=3.3

  crec=19.8

  k=19.6

Bem agora podemos chamar o modelo. Vamos rodá-lo e armazenar a saída em uma variável chamada result

#Definindo a área de drenagem
Area=4858
#Definindo o vetor de parâmetros
para=c(1795.2,3.3,19.8,19.6)

#Passando data frame para matrix
EPQ=as.matrix(EPQ)

#rodando o modelo e salvando a saída do processamento na variável result
result=smapM(EPQ,para,Area)

Bem, mas o que será que ficou gravando dentro de result? Vamos ver. Observe que temos 6 elementos e que para cada elemento temos um vetor contendo 57 elementos (exceto para vazões simuladas que temos 56 elementos, pois pulamos o primeiro no nosso loop for em função da inicialização das variáveis) . Como definimos anteriormente na nossa função, o elemento [[1]] é o vetor com as vazões simuladas, o elemento [[2]] o vetor com a série de nívels de Rsol, o elmento [[3]], o vetor com a série de níveis de rsub, o elemento [[4]], o vetor com a série de evapotranspiração real, o elemento [[5]], a série com os valores de escoamento superficial, e finalmente o elemento [[6]] a série com as vazões de base.

result
## [[1]]
##  [1]  26.06039  20.64469  22.35454  18.93134  27.34218  40.58950  74.92075
##  [8] 149.90246  61.37036  28.19990  56.22352  37.45133  24.52232  23.86259
## [15]  24.45223  34.05488  37.20387  37.13152  62.57816 160.67685 500.85353
## [22] 173.68091 260.05703  71.41083  56.41684  44.71798  43.63217  45.46013
## [29]  53.13275  66.76881  97.45833 165.88990 180.46778 132.38701 102.62306
## [36]  69.36586  72.23795  52.02511  58.68310  60.06270  48.93244  47.93335
## [43]  70.64577 162.13849 123.35341  63.61993 121.57784  88.16934  68.97231
## [50]  56.37912  49.41808  48.25445  59.75158  57.33550  74.77258 219.24330
## 
## [[2]]
##  [1]  718.0800  744.6420  723.3376  723.0845  684.0492  743.9677  852.8573
##  [8] 1011.1511 1179.4357 1124.0137 1031.3702 1041.3002 1012.5712  954.5803
## [15]  905.4582  861.2543  871.1182  885.4843  883.0540  971.8086 1175.3650
## [22] 1381.3693 1328.3722 1364.7739 1256.3371 1179.2520 1103.0504 1034.5509
## [29]  968.2178  941.1510  958.5080 1045.7612 1179.5783 1231.2262 1219.1506
## [36] 1175.4525 1112.4390 1084.8231 1025.7987  994.1570  966.8954  899.2576
## [43]  828.2908  891.9141 1132.2157 1146.2916 1077.1258 1135.0305 1119.6742
## [50] 1085.2367 1042.7161  981.1154  916.4901  911.5194  888.4225  946.1120
## [57] 1211.0295
## 
## [[3]]
##  [1] 311.6130 304.7133 298.2474 291.6484 284.8780 278.3786 274.6363
##  [8] 278.0950 299.1576 328.9423 345.4518 355.7911 365.4654 370.6243
## [15] 371.2015 368.7108 365.1049 362.3402 360.0899 360.2748 375.0042
## [22] 442.5262 513.8144 579.1913 634.0932 662.6021 676.8322 680.2051
## [29] 676.0277 667.6697 659.0641 655.0836 665.2286 692.5212 721.4086
## [36] 743.2456 754.9405 759.1091 758.0814 751.8209 743.0674 731.0323
## [43] 714.9734 698.7755 694.1886 708.7483 716.8144 722.6443 732.8465
## [50] 738.5591 739.1115 733.9733 723.3667 710.3017 696.9726 684.8371
## [57] 688.8714
## 
## [[4]]
##  [1]  0.00000 29.60325 27.85555 29.75320 38.09933 34.30660 47.03558
##  [8] 54.89713 65.07058 89.92905 82.12480 63.69036 48.83745 43.79959
## [15] 35.98604 37.01402 37.82474 46.88788 62.19503 63.12083 68.25810
## [22] 74.69194 90.67729 75.07054 68.22233 49.72444 43.15367 42.08922
## [29] 52.70632 49.51438 60.97220 64.05209 72.00020 87.61632 79.30817
## [36] 84.11151 66.62368 50.30846 42.31434 42.75885 50.87244 59.06903
## [43] 67.32280 66.24894 64.40189 86.37768 75.20450 72.70348 61.56959
## [50] 51.29629 39.01213 47.17977 54.27202 53.17974 67.69534 62.14455
## [57] 69.89659
## 
## [[5]]
##  [1]   0.00000000   3.28094637   0.58876102   1.73910207   0.11515539
##  [6]   4.90382543  12.30142428  31.01753588  71.49058916  22.82966048
## [11]   3.83708331  18.43472423   7.91267352   0.57707972   0.04066244
## [16]   0.33982617   5.62500423   7.45508085   7.51197760  21.36634791
## [21]  74.46810311 258.11947381  78.65021962 122.93506981  18.53509567
## [26]   8.51006576   1.18600814   0.10373017   0.97614670   5.27506182
## [31]  12.94769720  29.86125325  67.04670373  74.58630422  47.60824861
## [36]  30.49104850  11.72767013  12.87618867   1.78861432   5.42879262
## [41]   6.39320179   0.67171495   0.54901171  13.40294328  63.49763743
## [46]  42.65974749   9.81563086  40.91232641  22.62322512  11.87593947
## [51]   4.85980710   1.07208139   0.62065802   7.21346125   6.35941873
## [56]  16.26255695  94.89706825
## 
## [[6]]
##  [1]  0.000000 10.827501 10.587762 10.363090 10.133798  9.898549  9.672719
##  [8]  9.542687  9.662865 10.394719 11.429637 12.003286 12.362542 12.698692
## [15] 12.877947 12.897999 12.811457 12.686164 12.590099 12.511909 12.518335
## [22] 13.030132 15.376292 17.853318 20.124944 22.032602 23.023190 23.517639
## [29] 23.634834 23.489683 23.199273 22.900257 22.761949 23.114452 24.062777
## [36] 25.066518 25.825278 26.231636 26.376483 26.340774 26.123244 25.819089
## [43] 25.400910 24.842916 24.280095 24.120713 24.626613 24.906883 25.109454
## [50] 25.463947 25.662441 25.681633 25.503100 25.134555 24.680590 24.217450
## [57] 23.795783

EXPLORANDO GRAFICAMENTE OS RESULTADOS

Finalmente, vamos ver como ficou essa bagaceira de dados simulados. Será que as vazões simuladas ficaram próximas das vazões observadas? Antes de fazermos um gráfico bonito com o ggplot2, como estou curioso, vamos plotar rapidamente as vazões simuladas e as observadas para ver o que temos.

Como a série simulada ficou com um elemento a menos, pois pulamos um para inicialização, vamos comparar os valores simulados com os 56 primeiros elementos das vazões observaas, para mantermos compatíveis as comparações.Vamos salvar as vazões observadas na variável qobs e as vazões observadas na variável qsim.

qobs=EPQ[1:56,3]
qsim=result[[1]]

Agora vamos usar as funções padrões do R para plotar. Verifique a sintaxe para ver como plotamos linhas e inserimors legendas.

##Usando a função plot, passando o limite da variável y para melhor visualização, definindo título do gráfico, tipo de linha como ponto e tipo de ponto (pch=19), cor do ponto e títulos dos eixos x e y 
plot(qobs,ylim=c(0,500),main="Vazões simuladas X Vazões Observadas",type="p",pch=19,col="blue",xlab="Meses",ylab="Vazão (m³/s)")

#Passando um alinha azul do tipo tracejada (lty=2) em cima dos pontos observados
lines(qobs,col="blue",lty=2)

#Plotando a linha das vazões observadas
lines(qsim,col="red")

#Configurando a legenda
legend("topright", c("Observado", "Simulado"), col = c("blue","red"),
       lty = c(2,1), pch = c(19,NA))

Usando o ggplot 2 para gerar gráficos

Bem, alguns podem se dar por satisfeitos até aqui. Mas para quem quer ver o negócio ficar mais bacana, vamos em frente para melhorar a visualização dos nossos resultados e entender melhor o modelo usando o pacote ggplot2.

Faremos o seguinte, vamos criar uma estrutura chamada data frame (que nada mais é que uma variável do R que armazenará uma tabela podendo conter números e textos). Nesta tabela vamos inserir uma coluna com os meses que temos e mais duas colunas com os valores observados e simulados.O nosso número de meses será 56 no total. Abaixo o código comentado.

#Define o número total de meses
nmeses=56

#Cria uma sequencia com 56 meses iniciando no primeiro mês de maio de 1986
data=seq(as.Date("1933/5/1"), by = "month", length.out = nmeses)

#Vazões observadas
qobs=EPQ[1:56,3]

#Vazões simuladas
qsim=result[[1]]

#Juntando em uma so matriz as colunas das vazões observadas e simuladas com o comando cbind
vazoes=cbind(qobs,qsim)

#Forçando a matriz vazões a ser um data frame
vazoes=as.data.frame(vazoes)

#Inserindo mais uma coluna com as datas
dataFrame=cbind(data,vazoes)


#Dando nome as colunas do dataFrame
cols=c("Meses","Simulado","Observado")
colnames(dataFrame) = cols

Vamos ver como ficou esse negócio. Veja que agora temos os meses, as vazões simuladas e observadas em uma única estrutura de dados chamada dataFrame

#Define o número total de meses
nmeses=56

#Cria uma sequencia com 56 meses iniciando no primeiro mês de maio de 1986
data=seq(as.Date("1933/5/1"), by = "month", length.out = nmeses)

#Vazões observadas
qobs=EPQ[1:56,3]

#Vazões simuladas
qsim=result[[1]]

#Juntando em uma so matriz as colunas das vazões observadas e simuladas com o comando cbind
vazoes=cbind(qsim,qobs)

#Forçando a matriz vazões a ser um data frame
vazoes=as.data.frame(vazoes)

#Inserindo mais uma coluna com as datas
dataFrame=cbind(data,vazoes)


#Dando nome as colunas do dataFrame
cols=c("Meses","Simulado","Observado")
colnames(dataFrame) = cols

#Vendo como ficou nosso data Frame
head(dataFrame)
##        Meses Simulado Observado
## 1 1933-05-01 26.06039  88.70539
## 2 1933-06-01 20.64469  78.00904
## 3 1933-07-01 22.35454  66.49162
## 4 1933-08-01 18.93134  55.62024
## 5 1933-09-01 27.34218  61.90446
## 6 1933-10-01 40.58950  91.61406

Estamos quase prontos para gerar nosso gráfico com o ggplot2. Primeiramente vamos precisar agrupar nossos dados do data frame utilizando uma função do pacote reshape2, chamada melt. Antes de entender o que ela fez vamos executá-la com o comando a seguir.

dataFrameMelted <- reshape2::melt(dataFrame, id.var='Meses')

Bem, mas o que será que aconteceu? Veremos… (se chegou até aqui, é bom dar uma respirada, pois as vezes esse troca troca de formatação dos dados nos deixa meio tonto, mas no final vai compensar o sacrifício) .

Vemos a seguir que a nossa nova estrutura de dados dataFrameMelted criou duas colunas uma chamada variable e outra coluna chamada value. A coluna variable é como se fosse o tipo ou a classe da nossa variável, ou seja, se foi a vazão simulada é classificada como Simulado se for a vazão observada é classificada como Observado e a última coluna value mostra o respectivo valor da vazão.

dataFrameMelted <- reshape2::melt(dataFrame, id.var='Meses')

dataFrameMelted
##          Meses  variable     value
## 1   1933-05-01  Simulado  26.06039
## 2   1933-06-01  Simulado  20.64469
## 3   1933-07-01  Simulado  22.35454
## 4   1933-08-01  Simulado  18.93134
## 5   1933-09-01  Simulado  27.34218
## 6   1933-10-01  Simulado  40.58950
## 7   1933-11-01  Simulado  74.92075
## 8   1933-12-01  Simulado 149.90246
## 9   1934-01-01  Simulado  61.37036
## 10  1934-02-01  Simulado  28.19990
## 11  1934-03-01  Simulado  56.22352
## 12  1934-04-01  Simulado  37.45133
## 13  1934-05-01  Simulado  24.52232
## 14  1934-06-01  Simulado  23.86259
## 15  1934-07-01  Simulado  24.45223
## 16  1934-08-01  Simulado  34.05488
## 17  1934-09-01  Simulado  37.20387
## 18  1934-10-01  Simulado  37.13152
## 19  1934-11-01  Simulado  62.57816
## 20  1934-12-01  Simulado 160.67685
## 21  1935-01-01  Simulado 500.85353
## 22  1935-02-01  Simulado 173.68091
## 23  1935-03-01  Simulado 260.05703
## 24  1935-04-01  Simulado  71.41083
## 25  1935-05-01  Simulado  56.41684
## 26  1935-06-01  Simulado  44.71798
## 27  1935-07-01  Simulado  43.63217
## 28  1935-08-01  Simulado  45.46013
## 29  1935-09-01  Simulado  53.13275
## 30  1935-10-01  Simulado  66.76881
## 31  1935-11-01  Simulado  97.45833
## 32  1935-12-01  Simulado 165.88990
## 33  1936-01-01  Simulado 180.46778
## 34  1936-02-01  Simulado 132.38701
## 35  1936-03-01  Simulado 102.62306
## 36  1936-04-01  Simulado  69.36586
## 37  1936-05-01  Simulado  72.23795
## 38  1936-06-01  Simulado  52.02511
## 39  1936-07-01  Simulado  58.68310
## 40  1936-08-01  Simulado  60.06270
## 41  1936-09-01  Simulado  48.93244
## 42  1936-10-01  Simulado  47.93335
## 43  1936-11-01  Simulado  70.64577
## 44  1936-12-01  Simulado 162.13849
## 45  1937-01-01  Simulado 123.35341
## 46  1937-02-01  Simulado  63.61993
## 47  1937-03-01  Simulado 121.57784
## 48  1937-04-01  Simulado  88.16934
## 49  1937-05-01  Simulado  68.97231
## 50  1937-06-01  Simulado  56.37912
## 51  1937-07-01  Simulado  49.41808
## 52  1937-08-01  Simulado  48.25445
## 53  1937-09-01  Simulado  59.75158
## 54  1937-10-01  Simulado  57.33550
## 55  1937-11-01  Simulado  74.77258
## 56  1937-12-01  Simulado 219.24330
## 57  1933-05-01 Observado  88.70539
## 58  1933-06-01 Observado  78.00904
## 59  1933-07-01 Observado  66.49162
## 60  1933-08-01 Observado  55.62024
## 61  1933-09-01 Observado  61.90446
## 62  1933-10-01 Observado  91.61406
## 63  1933-11-01 Observado 139.12884
## 64  1933-12-01 Observado 274.52998
## 65  1934-01-01 Observado 111.01908
## 66  1934-02-01 Observado  70.83911
## 67  1934-03-01 Observado  73.16030
## 68  1934-04-01 Observado  66.26756
## 69  1934-05-01 Observado  49.60337
## 70  1934-06-01 Observado  42.01892
## 71  1934-07-01 Observado  39.41849
## 72  1934-08-01 Observado  42.17913
## 73  1934-09-01 Observado  47.39675
## 74  1934-10-01 Observado  41.78019
## 75  1934-11-01 Observado  68.43711
## 76  1934-12-01 Observado 174.76978
## 77  1935-01-01 Observado 392.44065
## 78  1935-02-01 Observado 217.04877
## 79  1935-03-01 Observado 274.65421
## 80  1935-04-01 Observado 134.03195
## 81  1935-05-01 Observado  90.89253
## 82  1935-06-01 Observado  72.89531
## 83  1935-07-01 Observado  61.82041
## 84  1935-08-01 Observado  54.57224
## 85  1935-09-01 Observado  52.70277
## 86  1935-10-01 Observado  54.23917
## 87  1935-11-01 Observado  83.92485
## 88  1935-12-01 Observado 121.87223
## 89  1936-01-01 Observado 189.88333
## 90  1936-02-01 Observado 131.57267
## 91  1936-03-01 Observado  91.08564
## 92  1936-04-01 Observado  65.29275
## 93  1936-05-01 Observado  64.05746
## 94  1936-06-01 Observado  49.93073
## 95  1936-07-01 Observado  45.33855
## 96  1936-08-01 Observado  40.68125
## 97  1936-09-01 Observado  32.12486
## 98  1936-10-01 Observado  29.29812
## 99  1936-11-01 Observado  39.25967
## 100 1936-12-01 Observado 108.41482
## 101 1937-01-01 Observado  94.81226
## 102 1937-02-01 Observado  55.01858
## 103 1937-03-01 Observado 100.70395
## 104 1937-04-01 Observado  69.10693
## 105 1937-05-01 Observado  50.25531
## 106 1937-06-01 Observado  44.95285
## 107 1937-07-01 Observado  35.66480
## 108 1937-08-01 Observado  30.24646
## 109 1937-09-01 Observado  33.93355
## 110 1937-10-01 Observado  31.39666
## 111 1937-11-01 Observado  41.68761
## 112 1937-12-01 Observado 197.04215

Uffa!!! Finalmente, vamos para o ggplot2!! Agora é um pouco mais chatinho para entender a fundo a sintaxe a primeira vez que vemos. Mas aos poucos dá para ir absorvendo. Deixei uns comentários no código para tentar dar uma ajudada.

##Chamando a função ggplot e passando como argumento o nosso data.FrameMelted, setando as variáveis x e y, e indicando que as cores, o tipo d elinha e o tamanho serão função do tipo ou classe da variável [fractor(variable)]. Salvei o resultado na variável "p" que é printada depois.

p<-ggplot(dataFrameMelted, 
       aes(x = Meses, y = value, 
           colour = factor(variable), linetype = factor(variable),
           size = factor(variable))) +
#Chamando o comando para traçar as linhas automaticamente
  geom_line()+
  
  #Chamando o comando para traçar os pontos e filtrando apenas para a classe de variáveis Observado
  geom_point(data = . %>% filter(variable == "Observado"),
             shape = 19, size = 1.5)+
  #Configurando manualmente as cores das linhas
  scale_color_manual("", values = c("red", "blue")) + 
  
    #Configurando manualmente os tipos de linha contínua (1) e tracejada(2)
  scale_linetype_manual("", values = c(1, 2))+
  
  #Configurando a espessura da linha manualmente
  scale_size_manual("", values=c(0.8,0.5))+
  
  #Defindo o título do eixo y e a posição legenda
  ylab("Vazões (m³/s)")+
  theme(legend.position="top")

print(p)

Usando o gganimate junto com o ggplot2 para gerar gráficos animados

Bem se você não achou que valeu a pena este esforço todo para usar o ggplot2, pois se deu por satisfeito com o visual do nosso gráfico convencional com o comando padrão plot, agora acho que você irá se convencer… Vamos ao clímax do nosso tutoral agora… Vamos animar este negócio…

Usaremos a funcão gganimate para animar nosso gráfico e ver a cada passo de tempo o valor simulado e o valor observado.

Antes, precisamos somar à variável p, que armazinou nosso ggplot anteriormente, a expressão transition_reveal passando como argumento a coluna Meses do nosso data frame.

E, por fim, chamar a função gganimate, passando o número de frames desejado e a resolução da imagem animada .gif que será gerada.

Vamos ver o que isso dá…

p <- p+transition_reveal(Meses) 

animate(p, fps = 12,height=600,width=900)

APÊNDICES: EXPLORANDO UM POUCO MAIS O MODELO COM OS GRÁFICOS ANIMADOS

Nossa, acho que este tutorial ficou muito grande :-/

Mas antes de finalizar acho que vale a pena usar os gráficos animados para dar uma olhada em outras informações que o modelo nos fornece, para não precisarmos fazer outro post só para isso. Se já está cansado deixe para dar uma olhada neste item depois.

Vamos usar os gráficos animados para ver, por exemplo, quanto da vazão total simulada o modelo estima que é escoamento de base e escoamento superficial. O gráfico a seguir mostra a vazão total e a vazão de base, a diferença é decorrente do escoamento superficial.

Vazão total e Escoamento de Base

#Inserindo as vazões de base no nosso data frame

#Vazão de bae
qbase=result[[6]]
qbase=qbase[1:56]
#Convertendo de mm para m³/s
qbase=qbase*Area/2630


vazoes=cbind(qsim,qbase,qobs)

#Forçando a matriz vazões a ser um data frame
vazoes=as.data.frame(vazoes)

#Inserindo mais uma coluna com as datas
dataFrame=0
dataFrame=cbind(data,vazoes)
dataFrame=as.data.frame(dataFrame)


#Dando nome as colunas do dataFrame
cols=c("Meses","Simulado","VazaoBase","Observado")
colnames(dataFrame) = cols

dataFrameMelted=0
dataFrameMelted <- reshape2::melt(dataFrame, id.var='Meses')


p<-ggplot(dataFrameMelted, 
       aes(x = Meses, y = value, 
           colour = factor(variable), linetype = factor(variable),
           size = factor(variable))) + 

  geom_area(data = . %>% filter(variable == "VazaoBase"),
            color = NA, fill = "lightblue", alpha = 0.5, show.legend = FALSE)+
  geom_line()+
  geom_point(data = . %>% filter(variable == "Observado"),
             shape = 19, size = 1.5)+

  scale_color_manual("", values = c("red", "lightblue", "black")) + 
  scale_linetype_manual("", values = c(4, 1, 1))+
  scale_size_manual("", values=c(0.5,1,1))+

  ylab("Vazões (m³/s)")+
  theme(legend.position="top")




p <- p+transition_reveal(Meses) 

animate(p, fps = 12)