Analises hidrologicos em R
II. Vazões mínimas
Neste módulo, vamos aprender a: 1. Ler dados de vazão com “url”. 2. Separar descarga em “baseflow” e “stormflow”. 3. Calcular vazões mínimas.
Será que a conversão de floresta para pastagem tem mudado as vazões mínimas?
Todos os arquivos apresentados aqui foram baixados do: http://www.snirh.gov.br/hidroweb/ O formato dos arquivos do hidroweb não podem ser lidos diretamente no R. Escrevi outro programa para reformatar os arquivos do hidroweb no formato csv.
1. Ler arquivos de vazão e fazer figuras
1.1. Ler arquivos diretamente da web.
Coloquei um arquivo csv no Google Dirve; pode lê-lo diretamente usando “read.csv” e “url”. “Qdf” é a variável que tem os dados. “df” significa “data.frame” (tabela).
Este arquivo é para o Rio Jamari em Ariquemes (15430000)
Qdf = read.csv(url("https://docs.google.com/spreadsheets/d/e/2PACX-1vTXp46RFu7HRh3uVgeWMWXIz_9zbotUlXRjd8JNCEoxsecbMRXIpfVJ9oLm3vHbVT06nZSxmKvLdaZu/pub?gid=586467714&single=true&output=csv"))Ou, você pode carregar o arquivo que você fez no Módulo 1:
pasta.dados = "C:/minha_pasta/Rcurso2021/dados/"
setwd(pasta.dados)
Qdf = read.csv("Q.Jaru.csv")E ver o que tem no arquivo nas primeiras 6 linhas
head(Qdf)## Date Q.cms Status
## 1 1977-12-01 NA NA
## 2 1977-12-02 NA NA
## 3 1977-12-03 NA NA
## 4 1977-12-04 NA NA
## 5 1977-12-05 NA NA
## 6 1977-12-06 NA NAe nas últimas 6 linhas:
tail(Qdf)## Date Q.cms Status
## 15181 2020-11-25 213.152 1
## 15182 2020-11-26 203.465 1
## 15183 2020-11-27 190.820 1
## 15184 2020-11-28 193.952 1
## 15185 2020-11-29 209.904 1
## 15186 2020-11-30 211.525 11.2. Formato para datas: as.Date
R trata datas como objetos específicos. Tem que formatar usando as.Date.
Qdf$Date = as.Date(Qdf$Date)
head(Qdf)## Date Q.cms Status
## 1 1977-12-01 NA NA
## 2 1977-12-02 NA NA
## 3 1977-12-03 NA NA
## 4 1977-12-04 NA NA
## 5 1977-12-05 NA NA
## 6 1977-12-06 NA NA1.3. Fazer uma figura de série temporal
plot(Qdf$Date,Qdf$Q.cms,type="l")
1.3.1 Faça os rótulos dos eixos.
plot(Qdf$Date,Qdf$Q.cms,type="l",xlab="Data",ylab="Vazao, metros cubicos por segundo")
1.3.2. Limitar os dados só para 2-3 anos:
datas.sub = as.Date(c("2006-08-01","2008-10-31"))
Qsub = Qdf[(Qdf$Date <= datas.sub[2]) & (Qdf$Date >= datas.sub[1]),]
plot(Qsub$Date,Qsub$Q.cms,type="l",xlab="Data",ylab="Vazao, metros cubicos por segundo")
2. Separação de vazão: pacote “EcohydRology”
R tem varios “libraries” ou grupos de programmas que tem que ser implementado. Use “install.packages” para baixar, e “library” para encargar o library.
Nos vamos usar “EcohydRology”, que tem varios funcoes para hidrologia. A descriccao de EcohydRology esta aqui: https://cran.r-project.org/web/packages/EcoHydRology/EcoHydRology.pdf
2.1. Instalar EcohydRology e hydroTSM
install.packages("EcohydRology")
install.packages("hydroTSM")2.2. Cargar EcohydRology
library(EcoHydRology)
library(hydroTSM)2.3. Ver a função “Baseflowseparation”
?BaseflowSeparation2.4. Faça a separação de vazão
2.4.1. A série não pode ter valores “NA”
q = Qsub$Q.cms
dates.nona = Qsub$Date[!is.na(q)]
q.nona = q[!is.na(q)]2.4.2. Implementar Baseflowseparation
bt eh “baseflow” qt eh “quickflow”
bsep = BaseflowSeparation(q.nona)
head(bsep)## bt qft
## 1 208.4367 21.24859
## 2 209.5058 37.16421
## 3 209.9444 33.28955
## 4 210.1806 29.63636
## 5 209.7294 30.08756
## 6 207.4277 28.99030tail(bsep)## bt qft
## 818 175.7915 47.21853
## 819 172.4361 43.98293
## 820 168.8497 60.82726
## 821 165.2367 44.66727
## 822 161.8761 41.58890
## 823 159.0167 31.803332.4.3. Coloca o baseflow na figura.
plot(Qsub$Date,Qsub$Q.cms,type="l",xlab="Data",ylab="Vazao, metros cubicos por segundo")
lines(dates.nona,bsep$bt,col="grey",lwd=2)
3. Vazões mínimas
3.1 Instalar libraries
Identificar vazões mínimas usando a library “hydrostats” e a função “low.spells”.
Só tem que instalar uma vez:
install.packages("hydrostats")E cargar:
library(hydrostats)## Warning: package 'hydrostats' was built under R version 4.1.13.2. Veja a função “low.spells”
?low.spellsquant=0.1 significa que vai achar os 10% valores menores
lows = low.spells(flow.ts = data.frame(Date=Qsub$Date,Q=Qsub$Q.cms),quant=0.1,plot=TRUE,ann.stats = TRUE)## Warning in max(x$lengths[which(x$values == 1)]): no non-missing arguments to
## max; returning -Inf## Warning in min(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## min; returning Inf## Warning in max(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## max; returning -Inf
# O que tem em “lows”?
lows## low.spell.threshold min.low.spell.duration avg.low.spell.duration
## 1 184.614 2 9
## med.low.spell.duration max.low.duration low.spell.freq avg.min.ann cv.min.ann
## 1 6 35 10 151.908 NA
## ann.min.timing ann.min.timing.sd ann.min.min.dur ann.min.avg.dur
## 1 264 0 Inf NaN
## ann.min.max.dur
## 1 -InfLimitar as análises para anos que têm dados para o ano inteiro (1971-2020):
Qdf.sub = Qdf[(Qdf$Date >= as.Date("1971-01-01")) & (Qdf$Date <= as.Date("2020-12-31")),]
Qdf.sub$Q = Qdf.sub$Q.cms # "low.spells" precisa que o nome do vazao ser "Q"3.3. Quantificar vazões mínimas para a série inteira:
3.4. Quantificar vazões mínimas para cada ano separado.
#Para isso, precisa de "plyr":
install.packages("plyr")library(plyr)## Warning: package 'plyr' was built under R version 4.1.1Faça uma coluna “ano”:
Qdf.sub$Ano = format(Qdf.sub$Date,"%Y")
# Qdf$Q = Qdf$Q.cms
anos = unique(Qdf.sub$Ano)Para cada ano, calcular a vazão mínima (10%), e colocar numa figura.
Qmins = data.frame(Ano=anos,Qmin=NA)
Qlen = data.frame(Ano=anos,Ndatas=NA)
for (i in 1:length(anos)){
ano = anos[i]
todosdatas = diy(ano)
Qsub = Qdf.sub[Qdf.sub$Ano==ano,]
Qlen[i,] = c(ano,length(!is.na(Qsub$Q.cms)))
Qinterp = approx(Qsub$Date, Qsub$Q, xout=as.Date(todosdatas))
Qinterp.df = data.frame(matrix(unlist(Qinterp),ncol=length(Qinterp),byrow=FALSE))
names(Qinterp.df) = c("Date","Q")
Qinterp.df$Date = as.Date(Qinterp.df$Date)
low.spells.data = low.spells(Qinterp.df, quant=0.1,plot=FALSE)
low.thresh = low.spells.data$low.spell.threshold
plot(Qinterp.df$Date,Qinterp.df$Q,type="l",main=ano,xlab="",ylab="Q cms")
lines(Qsub$Date,Qsub$Q,col="blue")
Qlow.values = Qinterp.df[Qinterp.df$Q<=low.thresh,]
points(Qlow.values$Date,Qlow.values$Q,col="red",pch=20)
Qmins$Qmin[i]= mean(Qlow.values$Q,na.rm=TRUE)
}## Warning in max(x$lengths[which(x$values == 1)]): no non-missing arguments to
## max; returning -Inf## Warning in min(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## min; returning Inf## Warning in max(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## max; returning -Inf
## Warning in max(x$lengths[which(x$values == 1)]): no non-missing arguments to
## max; returning -Inf## Warning in min(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## min; returning Inf## Warning in max(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## max; returning -Inf
## Warning in max(x$lengths[which(x$values == 1)]): no non-missing arguments to
## max; returning -Inf## Warning in min(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## min; returning Inf## Warning in max(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## max; returning -Inf
## Warning in max(x$lengths[which(x$values == 1)]): no non-missing arguments to
## max; returning -Inf## Warning in min(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## min; returning Inf## Warning in max(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## max; returning -Inf
## Warning in max(x$lengths[which(x$values == 1)]): no non-missing arguments to
## max; returning -Inf## Warning in min(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## min; returning Inf## Warning in max(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## max; returning -Inf
## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
## Warning in max(x$lengths[which(x$values == 1)]): no non-missing arguments to
## max; returning -Inf## Warning in min(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## min; returning Inf## Warning in max(ann.min.events.max.dur, na.rm = T): no non-missing arguments to
## max; returning -Inf
Qmins## Ano Qmin
## 1 1977 504.5680
## 2 1978 197.8170
## 3 1979 176.3579
## 4 1980 148.9213
## 5 1981 131.0626
## 6 1982 200.4470
## 7 1983 110.8060
## 8 1984 156.2814
## 9 1985 187.5778
## 10 1986 242.4044
## 11 1987 137.1843
## 12 1988 155.8036
## 13 1989 963.1379
## 14 1990 267.6119
## 15 1991 223.7045
## 16 1992 206.7573
## 17 1993 182.5533
## 18 1994 193.6643
## 19 1995 171.4238
## 20 1996 195.0835
## 21 1997 180.5380
## 22 1998 100.1150
## 23 1999 135.0359
## 24 2000 159.7374
## 25 2001 161.6262
## 26 2002 162.5200
## 27 2003 186.1578
## 28 2004 166.9223
## 29 2005 161.2029
## 30 2006 188.5087
## 31 2007 165.1982
## 32 2008 178.8712
## 33 2009 165.5627
## 34 2010 166.2533
## 35 2011 171.4454
## 36 2012 165.9003
## 37 2013 215.5777
## 38 2014 242.7976
## 39 2015 238.5201
## 40 2016 160.1874
## 41 2017 182.9277
## 42 2018 216.6706
## 43 2019 206.2848
## 44 2020 140.76493.5 Fazer um gráfico de vazões mínimas:
plot(Qmins$Ano,Qmins$Qmin,xlab="",ylab="Qmin, cms",type="l")
Incluir secas fortes identificado na literatura (Marengo et al., 2016)
anos.secas = c(1981,1983,1995,1998,2005,2010,2015)
Qmin.secas = merge(data.frame(Ano=anos.secas),Qmins,by.y="Ano")
plot(Qmins$Ano,Qmins$Qmin,xlab="",ylab="Qmin, cms",type="l")
points(Qmin.secas$Ano,Qmin.secas$Qmin,pch=20,col="red")
3.6 Identificar quais anos não têm dados completos.
Qmins$Qmin[as.numeric(Qlen$Ndatas)<365]=NA
plot(Qmins$Ano,Qmins$Qmin,xlab="",ylab="Qmin, cms",type="l")
points(Qmin.secas$Ano,Qmin.secas$Qmin,pch=20,col="red")