Analises hidrologicos em R
Todos os arquivos apresentado aqui foram baixado a http://www.snirh.gov.br/hidroweb/ O formato dos arquivos de hidroweb nao podem ser leito directamente em R. Escrevi outra programma para reformatar os arquivos no formato csv.
Como usar este documento
O codigo abaxio nas caixas cinzas pode ser copiado directamente em “Rstudio”, na janela para Scripts.
1. Ler arquivos de vazao e fazer figuras
1.1. Ler arquivos directamente do web.
Coloquei um arquivo csv em Google Dirve; pode le-lo directamente usando “read.csv” e “url”. “Qdf” eh a variavel que tem os dados. “df” significa “data.frame”.
Este arquivo eh 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"))E ver o que tem no arquivo nas primeiras 6 linhas
head(Qdf)## Date Q.cms
## 1 1970-07-01 NA
## 2 1970-07-02 NA
## 3 1970-07-03 NA
## 4 1970-07-04 NA
## 5 1970-07-05 NA
## 6 1970-07-06 NAe nas ultimas 6 linhas:
tail(Qdf)## Date Q.cms
## 16190 2016-01-26 96.73
## 16191 2016-01-27 71.53
## 16192 2016-01-28 71.53
## 16193 2016-01-29 78.20
## 16194 2016-01-30 88.28
## 16195 2016-01-31 77.231.2. Formato para datas: as.Date
R trata datas como objectos especificos. Tem que formatar usando as.Date.
Qdf$Date = as.Date(Qdf$Date)
head(Qdf)## Date Q.cms
## 1 1970-07-01 NA
## 2 1970-07-02 NA
## 3 1970-07-03 NA
## 4 1970-07-04 NA
## 5 1970-07-05 NA
## 6 1970-07-06 NA1.3. Fazer uma figura de serie temporal
plot(Qdf$Date,Qdf$Q.cms,type="l")
1.3.1 Faca os rotulos dos eixos.
plot(Qdf$Date,Qdf$Q.cms,type="l",xlab="Data",ylab="Vazao, metros cubicos por segundo")
1.3.2. Limitar os dados so 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. Separassao de vazao: 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
install.packages("EcohydRology")2.2. Cargar EcohydRology
library(EcoHydRology)2.3. Ver o funccao “Baseflowseparation”
?BaseflowSeparation2.4. Faca a separassao de vazao
2.4.1. A serie nao 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 16.33333 2.891378
## 2 13.02855 5.671453
## 3 12.60113 5.298868
## 4 12.19582 5.104182
## 5 11.79818 5.101818
## 6 11.41290 4.787101tail(bsep)## bt qft
## 755 5.557915 32.492085
## 756 6.014104 25.035896
## 757 6.422596 15.297404
## 758 6.780095 9.089905
## 759 7.089132 6.430868
## 760 7.357688 5.3323122.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. Vazoes minimas
3.1 Installar libraries
Identifica vazoes minimas usando library “hydrostats” e funccao “low.spells”.
So tem que installar uma vez:
install.packages("hydrostats")E cargar:
library(hydrostats)## Warning: package 'hydrostats' was built under R version 4.0.43.2. Veja a funccao “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)
# O que tem em “lows”?
lows## low.spell.threshold min.low.spell.duration avg.low.spell.duration
## 1 5.824 1 15.2
## med.low.spell.duration max.low.duration low.spell.freq avg.min.ann cv.min.ann
## 1 17 28 1.666667 NA NA
## ann.min.timing ann.min.timing.sd ann.min.min.dur ann.min.avg.dur
## 1 NA NA NA NA
## ann.min.max.dur
## 1 NALimitar as analyses para anos que tem dados para o ano enteiro (1971-2015):
Qdf.71.15 = Qdf[(Qdf$Date >= as.Date("1971-01-01")) & (Qdf$Date <= as.Date("2015-12-31")),]
Qdf.71.15$Q = Qdf.71.15$Q.cms # "low.spells" precisa que o nome do vazao ser "Q"3.3. Quantifica vazoes minimas para a serie enteira:
3.4. Quantifica vazoes minimas para cada ano separado.
#Para isso, precisa de "plyr":
install.packages("plyr")Faca uma coluna com “ano”:
Qdf.71.15$Ano = format(Qdf.71.15$Date,"%Y")
# Qdf$Q = Qdf$Q.cms
anos = unique(Qdf.71.15$Ano)Para cada ano, calcula a vazao minima (10%), e coloca numa figura.
Qmins = data.frame(Ano=anos,Qmin=NA)
for (i in 1:length(anos)){
ano = anos[i]
Qsub = Qdf.71.15[Qdf.71.15$Ano==ano,]
low.spells.data = low.spells(Qsub, quant=0.1,plot=FALSE)
low.thresh = low.spells.data$low.spell.threshold
plot(Qsub$Date,Qsub$Q,type="l",main=ano,xlab="",ylab="Q cms")
Qlow.values = Qsub[Qsub$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
Qmins## Ano Qmin
## 1 1971 2.9537838
## 2 1972 3.8010811
## 3 1973 6.2197368
## 4 1974 6.6537838
## 5 1975 5.6878947
## 6 1976 4.7160526
## 7 1977 11.0579487
## 8 1978 8.5255000
## 9 1979 12.2521053
## 10 1980 5.7982500
## 11 1981 4.4226316
## 12 1982 20.8869231
## 13 1983 8.5046154
## 14 1984 11.9294872
## 15 1985 7.8041026
## 16 1986 26.7924324
## 17 1987 7.3253846
## 18 1988 6.7692105
## 19 1989 12.8478378
## 20 1990 0.2840541
## 21 1991 0.2852778
## 22 1992 11.5991667
## 23 1993 9.9778947
## 24 1994 13.9240541
## 25 1995 11.1084615
## 26 1996 20.7394737
## 27 1997 20.8810811
## 28 1998 22.8382353
## 29 1999 11.0923913
## 30 2000 12.5989474
## 31 2001 27.5186486
## 32 2002 10.1585185
## 33 2003 17.8145161
## 34 2004 9.2474194
## 35 2005 4.1682857
## 36 2006 7.2470270
## 37 2007 4.3139394
## 38 2008 4.3744737
## 39 2009 6.0262500
## 40 2010 1.6565714
## 41 2011 1.4082857
## 42 2012 32.2697436
## 43 2013 2.0709091
## 44 2014 1.3762500
## 45 2015 0.12882353.5 Fazer um grafico de vazoes minimas:
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")