[R-br] Evapotranspiração Penman-Monteith
Yury Duarte
yurynepomuceno em gmail.com
Sexta Abril 29 14:12:32 BRT 2016
Boa tarde, Éder!!
Mais uma vez, lhe agradeço pelas contribuições. Foram de grande valia.
Como havia comentado, segue o script que programei para fazer o calculo
diário da ETP através dos métodos de Penman-Monteith (padronizado por
Allen, 1998) e através do método proposto por Priestly-Taylor.
Abraços!
rm(list = ls())
setwd("DIRETÓRIO")
a = read.table("ARQUIVO_CLIMATICO'); head(a)
a$Data = as.Date(a$Data)
a$DiaJuliano = as.numeric(format(a$Data, trim = T, '%j')); head(a)
altitude = 259.38
Patm = 101.3*((293-0.0065*altitude)/293)**5.26
gama = 0.665*(10**-3)*Patm
latitude = -7.53
corr = pi/180
decl = 23.45*sin(corr*((a$DiaJuliano-80)*360/365)); decl
hn = 1/corr*acos(-tan(corr*lat)*tan(corr*decl))
a$N = 2*hn/15
s = (1+0.033*cos(a$DiaJuliano*360/365))
t = corr*hn*sin(corr*lat)*sin(corr*decl)
u = cos(corr*lat)*cos(corr*decl)*sin(corr*hn)
a$Qo = 37.6*s*(t+u)
head(a)
x = 0.25
b = 0.50
n.est = (a$RS/a$Qo-x)*a$N/b; n.est
a$n<- ifelse((n.est)<0,1, n.est); head(a)
a$n.N = a$n/a$N; head(a)
r = 0.23
a$Tmed = (a$Tmax+a$Tmin)/2; head(a)
a$esMax = 0.6100*exp((17.3*a$Tmax)/(237.3+a$Tmax))
a$esMin = 0.6100*exp((17.3*a$Tmin)/(237.3+a$Tmin))
a$esMed = (a$esMax+a$esMin)/2
a$ea = (a$esMed*a$UR)/100
a$DPV = a$esMed-a$ea
a$Delta = (4098*a$esMed)/(a$Tmed+237.3)**2
a$BOC = a$RS*(1-r)
a$BOL =
-(0.903*(10**-9)*((a$Tmed+273)**4)*(0.34-0.14*(sqrt(a$ea)))*(0.1+(0.9*a$n.N)))
a$Rn = (a$BOC+a$BOL)
a$ETP.Penman =
((0.408*a$Delta*a$Rn)+(gama*(900/(a$Tmed+273))*a$Vel*a$DPV))/(a$Delta+(gama*(1+0.34*a$u2)))
#View(a)
a$W = ifelse(a$Tmed<16,(0.407+0.0145*a$Tmed),(0.483+0.01*a$Tmed))
a$ETP.Priestley = (1.26*a$W*a$Rn)/2.45
plot(a$ETP.Penman,a$ETP.Priestley)
write.table(a, "ETP.txt")
Yury Duarte
Engenheiro Agrônomo - ESALQ/USP
Em 28 de abril de 2016 14:51, Éder Comunello <comunello.eder em gmail.com>
escreveu:
> Yuri, boa tarde!
>
> Você deve editar as constantes conforme os dados disponíveis e fórmula
> utilizada.
>
> *Penman-Monteith diário utilizando "Rs":*
> # ET.PenmanMonteith(data, constants, ts="daily", solar="data", wind="yes",
> crop="short", ...)
> labels <- c("Elev", "lambda", "lat_rad", "Gsc", "z", "sigma", "G")
>
> *Penman-Monteith diário utilizando "n":*
> # ET.PenmanMonteith(data, constants, ts="daily", solar="sunshine hours",
> wind="yes", crop="short", ...)
> labels <- c("Elev", "lambda", "lat_rad", "Gsc", "z", "sigma", "G", "as",
> "bs")
>
> *Priestley-Taylor diário utilizando "Rs":*
> # ET.PriestleyTaylor(data, constants, ts="daily", solar="data",
> alpha=0.23, ...)
> labels <- c("Elev", "lambda", "lat_rad", "Gsc", "sigma", "G", "alphaPT")
>
> *Priestley-Taylor diário utilizando "n":*
> # ET.PriestleyTaylor(data, constants, ts="daily", solar="sunshine hours",
> alpha=0.23, ...)
> labels <- c("Elev", "lambda", "lat_rad", "Gsc", "sigma", "G", "alphaPT",
> "as", "bs")
>
> No caso de comparar resultados atente para diferenças de padrão das
> fórmulas para crop="short" (FAO-56) ou crop="tall" (ASCE-EWRI).
>
> Você pode avaliar diretamente as fórmulas empregadas no pacote! Se você
> digitar a função sem estar seguida de parênteses, terá acesso ao código
> utilizado e poderá comparar com teu script:
>
> > ET.PenmanMonteith
> # function (data, constants, ts = "daily", solar = "sunshine hours",
> # wind = "yes", crop = "short", ...)
> # {
> # if (is.null(data$Tmax) | is.null(data$Tmin)) {
> # stop("Required data missing for 'Tmax.daily' and 'Tmin.daily', or
> 'Temp.subdaily'")
> # }
> # if (is.null(data$RHmax) | is.null(data$RHmin)) {
> # stop("Required data missing for 'RHmax.daily' and 'RHmin.daily', or
> 'RH.subdaily'")
> # }
> # ...
>
>
> ================================================
> Éder Comunello
> Agronomist (UEM), MSc in Environ. Sciences (UEM)
> DSc in Agricultural Systems Engineering (USP/Esalq)
> Brazilian Agricultural Research Corporation (Embrapa)
> Dourados, MS, Brazil |<O>|
> ================================================
> GEO, -22.2752, -54.8182, 408m
> UTC-04:00 / DST: UTC-03:00
>
>
>
>
> Em 28 de abril de 2016 13:19, Yury Duarte <yurynepomuceno em gmail.com>
> escreveu:
>
>> Obrigado, Éder!!
>>
>> De fato, não consegui entender a necessidade de informar todas essas
>> constantes que vc abordou no código quando li o "help" do pacote.
>> Irei rodar seus comandos e acertarei as constantes segundo minhas
>> situações.
>> Estou finalizando a programação das fórmulas de PriestleyTaylor e
>> PenmanMonteith para comparar com as saídas das evapo fornecidas pelo
>> pacote. Assim que tiver os scripts prontos os enviarei nessa mesma conversa!
>>
>> Abs
>>
>> Yury Duarte
>> Engenheiro Agrônomo - ESALQ/USP
>>
>> Em 28 de abril de 2016 12:18, Éder Comunello <comunello.eder em gmail.com>
>> escreveu:
>>
>>> Yuri, bom dia!
>>>
>>> Esse pacote é um pouco melindroso quanto ao formato de entrada dos dados
>>> e muita coisa não está documentada.
>>>
>>> O objeto de entrada deve ser uma lista, iniciando com Date e J (data e
>>> dia juliano), sendo seguida pelas variáveis climáticas (na classe {zoo}).
>>>
>>> O script tá rodando a partir dos seus dados de exemplo, mas você tem que
>>> utilizar as constantes adequadas pro seu caso.
>>>
>>> ### <code r>
>>> rm(list=ls())
>>> require(Evapotranspiration)
>>>
>>> # URL <- "http://r-br.2285057.n4.nabble.com/attachment/4666053/0/a.txt"
>>>
>>> setwd("~/LAB/LEARN") ### Alterar!
>>> df <- read.table("a.txt", sep=";", head=T, as.is=T)
>>> df$Data <- as.Date(df$Data)
>>>
>>> climate <- lapply(as.list(df)[2:7], zoo, df$Data)
>>> J <- as.numeric(format(df$Data, "%j"))
>>> data <- c(list(Date.daily=df$Data, J=J), climate)
>>>
>>> names(data)
>>> # [1] "Date.daily" "J" "Rs" "Tmax" "Tmin"
>>> "RHmax" "u2" "RHmin"
>>> str(data)
>>>
>>> # Editar adequadamente as constantes utilizadas!!!
>>> # As constantes a serem utilizadas variam com a formulação escolhida...
>>> Ver ajuda...
>>> # ?ET.PenmanMonteith
>>> # ?ET.PristleyTaylor
>>> # pi/180*-23.45 # [1] -0.4092797
>>> myConst <- list(lambda = 2.45, sigma = 4.903e-09, Gsc = 0.082,
>>> lat = -23.45, lat_rad = -0.40928, as = 0.25,
>>> bs = 0.55, Elev = 480, z = 2, Roua = 1.2, Ca =
>>> 0.001013, G = 0,
>>> alphaA = 0.14, alphaPT = 1.26, ap = 2.4, fz =
>>> 28, b0 = 1,
>>> a_0 = 11.9, b_0 = -0.15, c_0 = -0.25, d_0 =
>>> -0.0107, e0 = 0.81917,
>>> e1 = -0.0040922, e2 = 1.0705, e3 = 0.065649,
>>> e4 = -0.0059684,
>>> e5 = -0.0005967, gammaps = 0.66, epsilonMo =
>>> 0.92, PA = 285.8,
>>> alphaMo = 17.27, betaMo = 237.3, sigmaMo =
>>> 5.67e-08, lambdaMo = 28.5,
>>> b1 = 14, b2 = 1.2)
>>>
>>> res1 <- ET.PenmanMonteith(data, myConst, ts="daily", solar="data",
>>> wind="yes", crop = "short")
>>> res2 <- ET.PriestleyTaylor(data, myConst, ts="daily", solar="data",
>>> alpha=.23)
>>>
>>> res <- cbind(PM=res1$ET.Daily, PT=res2$ET.Daily)
>>> head(res)
>>> fit <- lm(res$PM~res$PT-1); summary(fit)
>>> plot(res$PM~res$PT); abline(fit, col=2); abline(a=0, b=1, col=3, lty=2)
>>> # 1:1
>>> ### </code>
>>>
>>>
>>> ================================================
>>> Éder Comunello
>>> Agronomist (UEM), MSc in Environ. Sciences (UEM)
>>> DSc in Agricultural Systems Engineering (USP/Esalq)
>>> Brazilian Agricultural Research Corporation (Embrapa)
>>> Dourados, MS, Brazil |<O>|
>>> ================================================
>>> GEO, -22.2752, -54.8182, 408m
>>> UTC-04:00 / DST: UTC-03:00
>>>
>>>
>>>
>>>
>>> Em 27 de abril de 2016 18:09, Jônatan <jdtatsch em gmail.com> escreveu:
>>>
>>>> Qual erro? Qual os parâmetros usados no comando?
>>>> Exponha um exemplo reproduzível.
>>>>
>>>> 2016-04-27 15:02 GMT-03:00 Yury Duarte <yurynepomuceno em gmail.com>:
>>>>
>>>>> Boa tarde colegas programadores!
>>>>>
>>>>> Estou com algumas dificuldades em rodar o comando "ET.PenmanMonteith"
>>>>> do pacote "Evapotranspiration". Minha intenção é de calcular a
>>>>> evapotranspiração potencial para a escala diária. Segundo o pacote, as
>>>>> variáveis exigidas são:
>>>>> Tmax; Tmin; RHmax; RHmin; Rs e u2.
>>>>> Todas as variáveis exigidas pelo pacote estão compreendidas no meu
>>>>> arquivo de entrada (que segue em anexo), entretanto o comando segue dando
>>>>> erro.
>>>>> Caso alguém seja familiar com o pacote e/ou com o comando e puder me
>>>>> dar uma luz, seria de muita ajuda!
>>>>>
>>>>> Desde já agradeço pela colaboração de todos!
>>>>>
>>>>> Abraços
>>>>> Yury Duarte
>>>>> Engenheiro Agrônomo - ESALQ/USP
>>>>>
>>>>> _______________________________________________
>>>>> R-br mailing list
>>>>> R-br em listas.c3sl.ufpr.br
>>>>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>>>>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>>>>> código mínimo reproduzível.
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> ###############################################################
>>>> ## Jônatan Dupont Tatsch
>>>> ## Professor do Departamento de Física
>>>> ## Coordenador Substituto do Programa de Pós-Graduação em
>>>> Meteorologia (PPGMET)
>>>> ## Centro de Ciências Exatas e Naturais (CCNE)
>>>> ## Universidade Federal de Santa Maria - UFSM
>>>> ## Faixa de Camobi, Prédio 13 - Campus UFSM - Santa Maria, RS, Brasil
>>>> - 97105-900
>>>> ## Telefone: +55(55)33012083
>>>> ## www.ufsm.br/meteorologia
>>>> ###############################################################
>>>>
>>>> _______________________________________________
>>>> R-br mailing list
>>>> R-br em listas.c3sl.ufpr.br
>>>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>>>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>>>> código mínimo reproduzível.
>>>>
>>>
>>>
>>> _______________________________________________
>>> R-br mailing list
>>> R-br em listas.c3sl.ufpr.br
>>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>>> código mínimo reproduzível.
>>>
>>
>>
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160429/d7532428/attachment.html>
Mais detalhes sobre a lista de discussão R-br