In walmes/EACS: Estatística Aplicada à Ciência do Solo

source("config/setup.R")

Análise Exploratória

#-----------------------------------------------------------------------
# Carrega os pacotes necessários.

library(lattice)
library(latticeExtra)
library(nlme)
library(plyr)
library(wzRfun)
library(EACS)
library(captioner)

figs <- captioner(prefix = "Figura")
tbls <- captioner(prefix = "Tabela")

#-----------------------------------------------------------------------
# Análise exploratório dos dados.

# Pra facilitar o manuseio, vamos usar um nome curto.
cra <- teca_cra
str(cra)

cra$tens[cra$tens == 0] <- 0.1

cap <- figs(1, "Umidade do solo em função da tensão matricial\
                para 50 sítios cultivados com teca em 3\
                camadas do solo.")

xyplot(umid ~ tens | factor(loc), data = cra,
       groups = cam, type = c("o"), as.table = TRUE,
       strip = TRUE, layout = c(NA, 5),
       scales = list(x = list(log = 10)),
       xscale.components = xscale.components.log10ticks,
       auto.key = list(title = "Camada (cm)", cex.title = 1.1),
       ylab = expression("Umidade do solo" ~ (m^{3}~m^{-3})),
       xlab = expression(log[10] ~ "Tensão" ~ (kPa)))

#-----------------------------------------------------------------------
# Remove curvas atípicas.

del <- with(cra, {
    (loc == 4 & cam == levels(cam)[3]) |
        (loc == 27 & cam == levels(cam)[1]) |
        (loc == 37 & cam == levels(cam)[2]) |
        (loc == 47 & cam == levels(cam)[2]) |
        (loc == 47 & cam == levels(cam)[3])
})

cra <- droplevels(cra[!del, ])

Ajuste da Curva de Retenção de Água do Solo

O ajuste foi feito considerando a seguinte parametrização do modelo van Genuchten

$$ U(x) = U_r + \frac{U_s - U_r}{(1 + \exp{n(\alpha + x)})^{1 - 1/n}} $$

em que $U$ é umidade (m³ m^-3) do solo, $x$ é o log na base 10 da tensão matricial aplicada (kPa), $U_r$ é a umidade residual (assíntota inferior), $U_s$ é a umidade de satuação (assíntota superior), $\alpha$ e $n$ são parâmetros empíricos de forma da curva de retenção de água. Uma vez conhecido valores para as quatidades mencionadas, são obtidos

$$ \begin{align} S &= -n\cdot \frac{U_s - U_r}{(1 + 1/m)^{m + 1}} \newline I &= -\alpha - \log(m)/n \newline U_I &= U(x = I) \end{align} $$

em que $S$ é a taxa no ponto de inflexão, parâmetro que é tido como central para avaliação da qualidade física do solo, bem como $I$ que corresponde ao log da tensão no ponto de inflexão da curva de retenção de água do solo. A umidade correspondente a tensão no ponto de inflexão é representada por $U_I$.

O modelo foi ajustado aos dados usando o método de mínimos quadrados com algoritmo de Newton-Raphson para encontrar estimativas para os parâmetros.

#-----------------------------------------------------------------------
# Ajuste da CRA.

xyplot(umid ~ tens,
       # data = subset(cra, loc == 40 & cam == "[40, 80)"),
       data = cra,
       scales = list(x = list(log = 10)),
       xscale.components = xscale.components.log10ticks,
       ylab = expression("Umidade do solo" ~ (m^{3} ~ m^{-3})),
       xlab = expression(log[10] ~ "Tensão" ~ (kPa)))

# Logaritmo na base 10 da tensão matricial.
cra$ltens <- log10(cra$tens)

# Expressão do modelo van Genuchten.
model <- umid ~ Ur + (Us - Ur)/(1 + exp(n * (alp + ltens)))^(1 - 1/n)

# Valores iniciais para os parâmetros da curva.
start <- list(Ur = 0.3, Us = 0.6, alp = -0.5, n = 4)

n00 <- nls(model, data = cra, start = start)
coef(n00)

#-----------------------------------------------------------------------
# Ajustar para cada unidade experimental, 50 loc x 3 cam = 150 ue, se
# não tivesse sido removido algumas curvas.

cra$ue <- with(cra, interaction(loc, cam, drop = TRUE))
nlevels(cra$ue)

db <- groupedData(umid ~ ltens | ue,
                  data = cra, order.groups = FALSE)

n0 <- nlsList(model = model, data = db,
              start = as.list(coef(n00)))
c0 <- coef(n0)

pairs(c0)

# Alguma curva sem ajustar?
sum(!complete.cases(c0))

plot(augPred(n0),
     strip = FALSE,
     as.table = TRUE,
     ylab = expression("Umidade do solo" ~ (m^{3} ~ m^{-3})),
     xlab = expression(log[10] ~ "Tensão" ~ (kPa)))

#-----------------------------------------------------------------------
# Determinar os demais parâmetros da curva de água do solo.

params <- as.data.frame(
    do.call(rbind, strsplit(rownames(c0),
                            split = "\\.")))

names(params) <- c("loc", "cam")
params <- transform(params,
                    loc = as.integer(loc),
                    cam = factor(cam, levels = levels(cra$cam)))
params <- na.omit(cbind(params, c0))

params <- within(params, {
    m <- 1 - 1/n
    d <- Us - Ur
    S <- -d * n * (1 + 1/m)^(-m - 1)
    I <- -alp - log(m)/n
    Ui <- Ur + (Us - Ur)/(1 + exp(n * (alp + I)))^(1 - 1/n)
    cad <- Ui - Ur
    rm(d, m)
})

str(params)

# addmargins(xtabs(~loc + cam, data = params))
params <- arrange(params, loc, cam)

splom(params[, -(1:2)], type = c("p", "r"))

#-----------------------------------------------------------------------
# Salvando em um objeto no pacote.

teca_crapar <- params
str(teca_crapar)

Informações da sessão

cat(format(Sys.time(),
           format = "Atualizado em %d de %B de %Y.\n\n"))
sessionInfo()

walmes/EACS documentation built on Sept. 8, 2020, 12:50 p.m.