Determinação de variáveis derivadas
Roilan Hernández Valdés
Lendo os dados observados depois do procesamento de
all.vars <- readRDS("~/Dropbox/Dissertacao/ESTACOES/filled.est.MDC.rds")
info <- readRDS("~/Dropbox/Dissertacao/1_EstSP.MG_2001-2014_info.rds")
stacoes <- info$site[order(info$dist)]Calculando Pressão Parcial de Vapor de Água.
all.vars <- lapply(all.vars, function(i){
i <- mutate(.data = i,
es = rh /100 * 0.611 * exp( 17.27 * tair / (tair + 243.5)),
es.orig = rh.orig/100 * 0.611 * exp( 17.27 * tair.orig / (tair.orig + 243.5))
)
i
});names(all.vars) <- info$site ex.est <- 33
scatterPlot(all.vars[[ex.est]],
x = "tair",
y = "rh",
z = "es",
smooth = TRUE,
type = "daytime",
main = paste("UR vs TAR, niveis ES. Est",info$site[ex.est]),
ylab = "Umidade Relativa", xlab = "Temperatura do ar")
Determinando Radiação Potencial
all.vars <- lapply(info$site, function(i){
cat(i," .. ")
# Funçao do pacote REddyProc
b <- mutate(.data = all.vars[[i]],
potrad = as.vector(fCalcPotRadiation(DoY.V.n = as.numeric(format(all.vars[[i]][,"date"], "%j")),
Hour.V.n = as.numeric(format(all.vars[[i]][,"date"], "%H")),
Lat_deg.n = info$lat[which(info$site == i)],
Long_deg.n = info$lon[which(info$site == i)],
TimeZone_h.n = -3,
useSolartime.b = TRUE
))
)
b
});names(all.vars) <- info$site timePlot(selectByDate(all.vars[[ex.est]],year = 2009,month =7, day = 1:5),
c("rg","rg.orig","potrad"),group = TRUE,
ylab = "Radiação Global", xlab= "Hour",
lty = 1, lwd = 2,
main = paste("Radiação Global e Potencial, Julho. Est",
info$site[ex.est],"Lat",round(info$lat[ex.est],2) ),
date.breaks = 12*6, date.format = "%H")
timePlot(selectByDate(all.vars[[ex.est]],year = 2009,month =1, day = 1:5),
c("rg","rg.orig","potrad"),group = TRUE,
ylab = "Radiação Global", xlab= "Hour",
lty = 1, lwd = 2,
main = paste("Radiação Global e Potencial, Janeiro. Est",
info$site[ex.est],"Lat",round(info$lat[ex.est],2) ),
date.breaks = 12*6, date.format = "%H")
Determinando Radiação de Onda Longa incidente
Seleção do método para cálculo
Baseado em Marthews et al. (2011) escolho-se o método para cálculo da ROLi. A sugestão dada pelos autores é inviável pela falta dos dados de altura da base da nuvem nas estações automáticas. Decideu-se então a partir da Tabela 5 Marthews et al. (2011) a combinação de equações para obter \(\epsilon\) (emisividade de ar claro), \(\epsilon_e\) (emisividade efetiva do ar) e \(C\) (índice de cobertura de nuvens), que mostraram os menores RMS no periodo de maior sensibilidade. A combinação escolhida é:
Para horário diurno:
Angstrom (1915): \[\varepsilon = 0.83 - 0.18 * \left( 10^{-0.067*e} \right)\]
Jegede et al. (2006): \[ C =\begin{cases} 1.1 - \dfrac{H_{TOC}}{H_{CLR}} & 0 \le \dfrac{H_{TOC}}{H_{CLR}} \lt 0.9 \\ 2\left(1- \dfrac{H_{TOC}}{H_{CLR}} \right) & 0.9 \le \dfrac{H_{TOC}}{H_{CLR}} \le 1 \end{cases} \]
Angstrom (1915): \[\varepsilon_e = \varepsilon \left( 1 + 0.22* C \right)\]
Para horário noturno:
Brunt (1932) \[\varepsilon_a = 0.51 + 0.066 \sqrt e\]
\[L_{in} = \varepsilon_a \sigma T_a^4\]
Definindo períodos diurno e noturno
# NOTA:: Não consigo que funcionem todos os mutate ao mesmo tempo, "RSession abort"
# Definindo horários diurno e noturno
all.vars <- lapply(info$site, function(i){
b <- mutate(.data = all.vars[[i]],
daytime = ifelse(potrad > 0 , "day", "night")
)
b
});names(all.vars) <- info$siteCorreção a valores de Rg a partir da Radiação Potencial
# Corrigindo os valores de RG superiores a Radiação Potencial
all.vars <- lapply(info$site, function(i){
b <- mutate(.data = all.vars[[i]],
rg = ifelse(rg > potrad, potrad, rg)
)
b
});names(all.vars) <- info$site
# Corrigindo os valores de RG superiores a Radiação Potencial
all.vars <- lapply(info$site, function(i){
b <- mutate(.data = all.vars[[i]],
rg.orig = ifelse(rg.orig > potrad, potrad, rg.orig)
)
b
});names(all.vars) <- info$siteÍndice de “Clearness”
# Clearness indice
tmp.zlwd <- lapply(info$site, function(i){
ind.k <- all.vars[[i]][,"rg"] / all.vars[[i]][,"potrad"]
ind.k.orig <- all.vars[[i]][,"rg.orig"] / all.vars[[i]][,"potrad"]
return(data.frame(date = all.vars[[i]][,"date"],
daytime = all.vars[[i]][,"daytime"],
ind.k = ind.k,
ind.k.orig = ind.k.orig))
});names(tmp.zlwd) <- info$site timePlot(selectByDate(tmp.zlwd[[ex.est]], year = 2012, month = 5:7),
c("ind.k","ind.k.orig"),group = TRUE,
lty = c(1,1), col = c("red","green"),
date.breaks = 12, date.format = "%Y\n%m",
main = "Razão entre o Rg e Radiação Potencial"
)
Índice de cobertura de nuvens
# Índice de cobertura de nuvens
tmp.zlwd <- lapply(info$site, function(i){
b <- mutate(.data = tmp.zlwd[[i]],
ind.c = ifelse(ind.k <0.9, 1.1 - ind.k, 2*(1-ind.k))
)
b
});names(tmp.zlwd) <- info$site
# NOTA:: RSession aborta se uso os dois mutate ao mesmo tempo
tmp.zlwd <- lapply(info$site, function(i){
b <- mutate(.data = tmp.zlwd[[i]],
ind.c.orig = ifelse(ind.k.orig <0.9, 1.1 - ind.k.orig, 2*(1-ind.k.orig))
)
b
});names(tmp.zlwd) <- info$site timePlot(selectByDate(tmp.zlwd[[3]], year = 2012, month = 5:7),
c("ind.c","ind.c.orig"),
lty = c(1,1), col = c("red","green"),
group = TRUE,
date.breaks = 12, date.format = "%Y\n%m",
main = "Indice de cobertura de nuvens"
)
Emisividade e emisividade efetiva
# Emissividade Angstrom (1915)
tmp.zlwd <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = tmp.zlwd[[i]],
emiss_CAN = 0.83 - 0.18 * (10^(-0.067*all.vars[[i]][,"es"]) )
)
b
}); names(tmp.zlwd) <- info$site
tmp.zlwd <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = tmp.zlwd[[i]],
emiss_CAN.orig = 0.83 - 0.18 * (10^(-0.067*all.vars[[i]][,"es.orig"]) )
)
b
}); names(tmp.zlwd) <- info$site
# Emissividade Brunt (1932)
tmp.zlwd <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = tmp.zlwd[[i]],
emiss_CBR = 0.51 + 0.066*sqrt(all.vars[[i]][,"es"])
)
b
}); names(tmp.zlwd) <- info$site
tmp.zlwd <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = tmp.zlwd[[i]],
emiss_CBR.orig = 0.51 + 0.066*sqrt(all.vars[[i]][,"es.orig"])
)
b
}); names(tmp.zlwd) <- info$site scatterPlot(tmp.zlwd[[ex.est]],
x = "emiss_CAN",
y = "emiss_CBR",
main = "Emissividade Angstrom (1915) vs Brunt (1932) Preenchido")
scatterPlot(tmp.zlwd[[ex.est]],
x = "emiss_CAN.orig",
y = "emiss_CBR.orig",
main = "Emissividade Angstrom (1915) vs Brunt (1932) Original")
# Emissividade efetiva para horário diurno Angstrom (1915)
tmp.zlwd <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = tmp.zlwd[[i]],
emiss_AAN = emiss_CAN*(1+0.22*ind.c)
)
b
}); names(tmp.zlwd) <- info$site
tmp.zlwd <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = tmp.zlwd[[i]],
emiss_AAN.orig = emiss_CAN.orig*(1+0.22*ind.c.orig)
)
b
}); names(tmp.zlwd) <- info$site timePlot(selectByDate(tmp.zlwd[[ex.est]], year = 2008:2011,month = 1:7),
c("emiss_AAN","emiss_AAN.orig"),
group = TRUE, stack = TRUE,
lty = c(1,1), col = c("red", "green"),
main = "Emissividade Efetiva; Angstrom (1915) ")
scatterPlot(tmp.zlwd[[88]],
x = "emiss_CAN",
y = "emiss_AAN",
z = "ind.c",
smooth = TRUE)
Radiação de Onda Longa incidente
#Constante de Stefan-Boltzman [W m-2 K-4]
sigma <- 5.67*10^-8
# Cáculo da Radiação de Onda Longa incidente
# a partir da formulação diferenciada para periodos diurno e noturno
# Dia:
# ---- emissividade: eps(es) CAN Angstrom (1915)
# ---- emissividade efetiva: eps_e(eps,c) AAN Angstrom (1915)
# ---- cobertura de nuvens: c(Rg, RgPot) Jegede et al. (2006)
# Noite:
# ---- emissividade: eps_e(es) CBR Brunt (1932)
# ^^^^^
# zlwd = eps_e * sigma * Tair^4
# .........................................................................................
# Radiação de Onda Longa com os dados observados SEM preencimento de falhas
tmp.zlwd <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = tmp.zlwd[[i]],
zlwd = ifelse(daytime == "day",
emiss_AAN * sigma * (all.vars[[i]][,"tair"]+273.15)^4 ,
emiss_CBR * sigma * (all.vars[[i]][,"tair"]+273.15)^4 )
)
b
}); names(tmp.zlwd) <- info$site
# Radiação de Onda Longa com os dados observados COM preencimento de falhas
tmp.zlwd <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = tmp.zlwd[[i]],
zlwd.orig = ifelse(daytime == "day",
emiss_AAN.orig * sigma * (all.vars[[i]][,"tair.orig"]+273.15)^4 ,
emiss_CBR.orig * sigma * (all.vars[[i]][,"tair.orig"]+273.15)^4 )
)
b
}); names(tmp.zlwd) <- info$site timePlot(tmp.zlwd[[3]], c("zlwd","zlwd.orig"))
scatterPlot(tmp.zlwd[[ex.est]],
x = "ind.c",
y = "zlwd",
z = "ind.k",
type = "daytime", smooth = TRUE)
monthly.zlwd <- timeVariation(tmp.zlwd[["X001"]],
c("zlwd.orig","zlwd" ),
type = "month") update(monthly.zlwd$plot$hour, layout = c(6,2))
### TODO Pasar a função em script
require(Hmisc)
var <- "zlwd"
probs <- c(0.001,0.999)
ex.est <- 88
hours <- ifelse((combhur <- 0:23) < 10,paste("0",combhur,sep = ""), combhur )
by.hours <- sapply(hours, FUN = function(i){ #i <- "04"
x <- tmp.zlwd[[ex.est]][which(format(tmp.zlwd[[ex.est]][,"date"], "%H") == i),var]
x
})
by.hours <- as.data.frame(by.hours)
n.obs.by.hh <- max(sapply(X = by.hours,length))
by.hours <- lapply(by.hours, function(i){c(i,rep(NA, times = n.obs.by.hh - length(i)))})
by.hours <- as.data.frame(by.hours)
up.lim <- apply(by.hours, 2,quantile, probs = probs[2], na.rm = TRUE)
down.lim <- apply(by.hours, 2,quantile, probs = probs[1], na.rm = TRUE)
lty.f <- c(5,4)
lwd.f <- c(1,1)
boxplot(by.hours,
pch = 20 , cex = 0.5,col = "lightblue",
axes = FALSE,
ylab = "Valor" ,
xlab = "Horas",
main = paste("Estação", info$site[ex.est], var))
box(which = "plot")
axis(1, at = 1:24,labels = hours, tick = TRUE, padj = 0)
axis(2)
lines(x = down.lim, lwd = lwd.f , lty = lty.f[1], col = "red")
lines(x = up.lim, lwd = lwd.f, lty = lty.f[2], col = "purple")
legend("topright",
xpd = TRUE,
legend = c(paste("Percentil",probs[1]*100),paste("Percentil",probs[2]*100)),
bg = "transparent",
bty = "n",
lty = lty.f,
lwd = lwd.f,
inset = c( 0.0,-0.1)) 
Operações finais
# Escrevendo variável Radiação de Onda Longa no data frame geral
all.vars <- lapply(info$site, function(i){
cat(".",i,".")
b <- mutate(.data = all.vars[[i]],
zlwd = tmp.zlwd[[i]][,"zlwd"],
zlwd.orig = tmp.zlwd[[i]][,"zlwd.orig"])
b
}); names(all.vars) <- info$site # Salvando todas as variáveis a um arquivo
saveRDS(all.vars, "~/Dropbox/Dissertacao/EstSP.MG_2001-2014_FINAL.rds")Avaliações
Componentes de radiação.
yearly.all <- timeVariation(all.vars[["X001"]] ,
c("zlwd","zlwd.orig","rg","rg.orig","potrad" ),
type = "year")Ciclo médio diário por ano da Radiação Potencial, Radiação Global e Radiação de Onda Longa Incidente.
update(yearly.all$plot$hour, layout = c(5,3))
Ciclo médio mensal por ano da Radiação Potencial, Radiação Global e Radiação de Onda Longa Incidente.
update(yearly.all$plot$month, layout = c(5,3))