Preparação de dados horários de NEE observados sobre Cerrado
Jônatan Tatsch (Departamento de Física - UFSM)
31-05-2015
- 1 Introdução
- 2 Dados
- 3 Software e sistema operacional
- 4 Pacotes necessários e scripts
- 5 Preparação de dados
- 5.1 Importação dos dados
- 5.2 Formatação de colunas, nome das variáveis e ajuste de unidades
- 5.3 Adição de coluna com umidade do solo superficial
- 5.4 Adição de coluna DateTime
- 5.5 Adicionando coluna com Radiaçao solar Global incidente
- 5.6 Verificação da sequência cronológica dos dados
- 5.7 Conversão dos dados para classe de dados do REddyProc
- 5.8 Remoção ou Substituição de dados inconsistentes
- 6 Salvando dados de entrada para REddyProc
- 7 Variação temporal das variáveis
1 Introdução
Neste documento descreve-se a preparação dos dados horários de NEE para o preenchimento de falhas do NEE noturno (Respiração do ecossistema) em condições de turbulência suficiente a partir da relação com variáveis ambientais (temperatura do ar, déficit da pressão de vapor d’água, umidade do solo) e biofísicas (LAI) para o ecossistema Cerrado.
A preparação consistiu nas seguintes etapas de adaptação do conjunto de dados à estrutura e formato de dados requerido pelo pacote REddyProc:
ajuste de nomes e unidades das variáveis
inclusão dos dados de umidade do solo
adição de variável indicativa do tempo (
DateTime) no formatoPOSIXtque inclui informação condensada de ano, mês, dia, hora e minuto.inclusão de dados de radiação solar global incidente (
Rg) via regressão comPARverificação de sequência cronológica dos dados (dados igualmente espaçados temporalmente)
teste de plausibilidade dos valores das variáveis (via
REddyProc)visualização das séries temporais das variáveis
2 Dados
Os dados foram fornecidos em formato Excel (arquivo Jonatan_PdG_NEE_2009_2014_30min.xlsx) pelo Pesquisador Dr. Osvaldo Cabral da Embrapa Meio Ambiente.
Fluxos médios de 30 min de NEE foram medidos por um sistema de covariância dos vórtices turbulentos sobre o dossel do cerrado no sudeste do Brasil (Cabral et. al. 2015, In prep.).
Esses dados foram filtrados para u_star > 0.2 m s-1 e corrigidos pelo termo de armazenamento de CO2 no dossel. Nos períodos filtrados o NEE foi substituído por -9999.999. Esse pré-processamento dos dados de NEE foi realizado por Cabral et. al. 2015 (In Prep.).
As falhas de umidade do solo foram simuladas pelo modelo Hydrus 1D (Simunek et al. 2008). Os valores de 30 min foram repetidos dos valores diários (48 valores por dia, 30 min).
3 Software e sistema operacional
Os procedimento aqui descritos foram realizados no software R. Os códigos utilizados e o procedimentos necessários para aplicação daquelas técnicas são descritos. Os scripts foram criados e testados no Ambiente Integrado de Desenvolvimento (IDE) RStudio. O sistema operacional é Linux Ubuntu 14.04.2 LTS.
4 Pacotes necessários e scripts
Limpeza do espaço de trabalho e diretório de trabalho.
## clean any previous data
rm(list = ls())
## work directory
#getwd()
## diretórios
#dir()
options(digits = 6)Carregando pacotes necessários.
pcks <- c("plyr", "dplyr","reshape2","lubridate","stringr",
"openair","psych","REddyProc","ggplot2","scales",
"lattice", "latticeExtra")
invisible(sapply(pcks, require, character.only = TRUE, quietly = TRUE))Carregando scripts desenvolvidos.
source("R/change_hour_string.R")5 Preparação de dados
5.1 Importação dos dados
Os dados foram lidos no R extraindo primeiro a matriz com os valores e depois os nome das variáveis.
nee_d <- read.csv(file = "data/Jonatan_PdG_NEE_2009_2014_30min.csv",
head = TRUE,
na.strings = "-9999.999")
names(nee_d) <- tolower(names(nee_d))
names(nee_d) <- gsub("umol", "", names(nee_d))
names(nee_d) <- gsub("_kpa", "", names(nee_d))
names(nee_d) <- gsub("_c", "", names(nee_d))
names(nee_d) <- gsub("dfct", "vpd", names(nee_d))
names(nee_d) <- gsub("par", "PAR", names(nee_d))
head(nee_d)| day | year | doy | hour | nee | PAR | ustar | vpd | tar | lai | swc1m | swc2m |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 189.021 | 2009 | 189 | 50 | 2.11495 | 0 | 0.580 | 0.422681 | 17.370 | 1.5226 | 46.7653 | 123.205 |
| 189.042 | 2009 | 189 | 100 | 4.38557 | 0 | 0.645 | 0.380730 | 16.943 | 1.5226 | 46.7653 | 123.205 |
| 189.062 | 2009 | 189 | 150 | 2.11495 | 0 | 0.618 | 0.381223 | 16.870 | 1.5226 | 46.7653 | 123.205 |
| 189.083 | 2009 | 189 | 200 | 2.40042 | 0 | 0.533 | 0.367781 | 16.700 | 1.5226 | 46.7653 | 123.205 |
| 189.104 | 2009 | 189 | 250 | 1.73946 | 0 | 0.520 | 0.371807 | 16.653 | 1.5226 | 46.7653 | 123.205 |
| 189.125 | 2009 | 189 | 300 | 3.13241 | 0 | 0.504 | 0.364160 | 16.537 | 1.5226 | 46.7653 | 123.205 |
5.2 Formatação de colunas, nome das variáveis e ajuste de unidades
Ordenando colunas e ajustando nomes das variáveis conforme arquivo ASCII de exemplo do pacote REddyProc:
variáveis: Year DoY Hour (NEE LE H Rg Tair Tsoil rH VPD Ustar)
unidades: (umolm-2s-1 Wm-2 Wm-2 Wm-2 degC degC % hPa ms-1)
Uma descrição sobre o formato dos dados de entrada estão disponíveis aqui.
nee_d <- subset(nee_d, sel = c("year", "doy", "hour", "nee", "PAR", "tar",
"vpd", "ustar", "lai", "swc1m", "swc2m"))
names(nee_d) <- c("Year","DoY", "Hour", "NEE", "PAR", "Tair", "VPD", "Ustar",
"lai", "swc1m", "swc2m")
str(nee_d)'data.frame': 87360 obs. of 11 variables:
$ Year : int 2009 2009 2009 2009 2009 2009 2009 2009 2009 2009 ...
$ DoY : int 189 189 189 189 189 189 189 189 189 189 ...
$ Hour : int 50 100 150 200 250 300 350 400 450 500 ...
$ NEE : num 2.11 4.39 2.11 2.4 1.74 ...
$ PAR : num 0 0 0 0 0 0 0 0 0 0 ...
$ Tair : num 17.4 16.9 16.9 16.7 16.7 ...
$ VPD : num 0.423 0.381 0.381 0.368 0.372 ...
$ Ustar: num 0.58 0.645 0.618 0.533 0.52 0.504 0.498 0.512 0.591 0.601 ...
$ lai : num 1.52 1.52 1.52 1.52 1.52 ...
$ swc1m: num 46.8 46.8 46.8 46.8 46.8 ...
$ swc2m: num 123 123 123 123 123 ...# convertendo VPD de kPa par hPa
nee_d$VPD <- nee_d$VPD*10
head(nee_d)| Year | DoY | Hour | NEE | PAR | Tair | VPD | Ustar | lai | swc1m | swc2m |
|---|---|---|---|---|---|---|---|---|---|---|
| 2009 | 189 | 50 | 2.11495 | 0 | 17.370 | 4.22681 | 0.580 | 1.5226 | 46.7653 | 123.205 |
| 2009 | 189 | 100 | 4.38557 | 0 | 16.943 | 3.80730 | 0.645 | 1.5226 | 46.7653 | 123.205 |
| 2009 | 189 | 150 | 2.11495 | 0 | 16.870 | 3.81223 | 0.618 | 1.5226 | 46.7653 | 123.205 |
| 2009 | 189 | 200 | 2.40042 | 0 | 16.700 | 3.67781 | 0.533 | 1.5226 | 46.7653 | 123.205 |
| 2009 | 189 | 250 | 1.73946 | 0 | 16.653 | 3.71807 | 0.520 | 1.5226 | 46.7653 | 123.205 |
| 2009 | 189 | 300 | 3.13241 | 0 | 16.537 | 3.64160 | 0.504 | 1.5226 | 46.7653 | 123.205 |
5.3 Adição de coluna com umidade do solo superficial
Dados de armazenamento de água na camada superficial do solo (0-21 cm) na escala diária foram fornecidos em um arquivo separado, gerado previamente quando haviam somente os dados diários.
No trecho de código abaixo vamos combinar os dados diários (com swctop) com os de 30 min (nee_d).
Importando dados de umidade do solo:
swc_top <- read.csv2("data/Jonatan_Umidade_Hydrus1.csv", head = TRUE)
names(swc_top)[names(swc_top) == "X10.5cm"] <- "swc_11cm"
names(swc_top)[names(swc_top) == "X21cm"] <- "swc_21cm"Corrigindo coluna Ano que está com valores do ano anterior nos 5 primeiros dias do ano atual equivoco.
# onde está o erro
subset(swc_top, DOY %in% 1:7)| Ano | DOY | swc_11cm | swc_21cm | |
|---|---|---|---|---|
| 184 | 2009 | 1 | 0.1113 | 0.1206 |
| 185 | 2009 | 2 | 0.1018 | 0.1112 |
| 186 | 2009 | 3 | 0.0942 | 0.1021 |
| 187 | 2009 | 4 | 0.0870 | 0.0950 |
| 188 | 2009 | 5 | 0.0837 | 0.0887 |
| 189 | 2010 | 6 | 0.0860 | 0.0849 |
| 190 | 2010 | 7 | 0.1362 | 0.0839 |
| 549 | 2010 | 1 | 0.1400 | 0.1032 |
| 550 | 2010 | 2 | 0.1983 | 0.1982 |
| 551 | 2010 | 3 | 0.1926 | 0.1926 |
| 552 | 2010 | 4 | 0.1381 | 0.1457 |
| 553 | 2010 | 5 | 0.1173 | 0.1278 |
| 554 | 2010 | 6 | 0.1183 | 0.1176 |
| 555 | 2011 | 7 | 0.1937 | 0.1937 |
| 914 | 2011 | 1 | 0.1978 | 0.1978 |
| 915 | 2011 | 2 | 0.1281 | 0.1406 |
| 916 | 2011 | 3 | 0.1104 | 0.1208 |
| 917 | 2011 | 4 | 0.0992 | 0.1092 |
| 918 | 2011 | 5 | 0.1496 | 0.1125 |
| 919 | 2011 | 6 | 0.1495 | 0.1468 |
| 920 | 2012 | 7 | 0.1232 | 0.1329 |
| 1280 | 2012 | 1 | 0.1101 | 0.1191 |
| 1281 | 2012 | 2 | 0.1005 | 0.1098 |
| 1282 | 2012 | 3 | 0.1965 | 0.1957 |
| 1283 | 2012 | 4 | 0.1279 | 0.1404 |
| 1284 | 2012 | 5 | 0.1104 | 0.1207 |
| 1285 | 2013 | 6 | 0.1103 | 0.1110 |
| 1286 | 2013 | 7 | 0.1319 | 0.1135 |
| 1645 | 2013 | 1 | 0.1508 | 0.1496 |
| 1646 | 2013 | 2 | 0.1220 | 0.1325 |
| 1647 | 2013 | 3 | 0.1065 | 0.1164 |
| 1648 | 2013 | 4 | 0.0961 | 0.1051 |
| 1649 | 2013 | 5 | 0.0880 | 0.0967 |
| 1650 | 2014 | 6 | 0.0812 | 0.0892 |
| 1651 | 2014 | 7 | 0.0749 | 0.0827 |
# correções
swc_top[184:188,"Ano"] <- 2010
swc_top[549:554,"Ano"] <- 2011
swc_top[914:919,"Ano"] <- 2012
swc_top[1280:1284,"Ano"] <- 2013
swc_top[1645:1649,"Ano"] <- 2014
# verificando correção
subset(swc_top, DOY %in% 1:6)| Ano | DOY | swc_11cm | swc_21cm | |
|---|---|---|---|---|
| 184 | 2010 | 1 | 0.1113 | 0.1206 |
| 185 | 2010 | 2 | 0.1018 | 0.1112 |
| 186 | 2010 | 3 | 0.0942 | 0.1021 |
| 187 | 2010 | 4 | 0.0870 | 0.0950 |
| 188 | 2010 | 5 | 0.0837 | 0.0887 |
| 189 | 2010 | 6 | 0.0860 | 0.0849 |
| 549 | 2011 | 1 | 0.1400 | 0.1032 |
| 550 | 2011 | 2 | 0.1983 | 0.1982 |
| 551 | 2011 | 3 | 0.1926 | 0.1926 |
| 552 | 2011 | 4 | 0.1381 | 0.1457 |
| 553 | 2011 | 5 | 0.1173 | 0.1278 |
| 554 | 2011 | 6 | 0.1183 | 0.1176 |
| 914 | 2012 | 1 | 0.1978 | 0.1978 |
| 915 | 2012 | 2 | 0.1281 | 0.1406 |
| 916 | 2012 | 3 | 0.1104 | 0.1208 |
| 917 | 2012 | 4 | 0.0992 | 0.1092 |
| 918 | 2012 | 5 | 0.1496 | 0.1125 |
| 919 | 2012 | 6 | 0.1495 | 0.1468 |
| 1280 | 2013 | 1 | 0.1101 | 0.1191 |
| 1281 | 2013 | 2 | 0.1005 | 0.1098 |
| 1282 | 2013 | 3 | 0.1965 | 0.1957 |
| 1283 | 2013 | 4 | 0.1279 | 0.1404 |
| 1284 | 2013 | 5 | 0.1104 | 0.1207 |
| 1285 | 2013 | 6 | 0.1103 | 0.1110 |
| 1645 | 2014 | 1 | 0.1508 | 0.1496 |
| 1646 | 2014 | 2 | 0.1220 | 0.1325 |
| 1647 | 2014 | 3 | 0.1065 | 0.1164 |
| 1648 | 2014 | 4 | 0.0961 | 0.1051 |
| 1649 | 2014 | 5 | 0.0880 | 0.0967 |
| 1650 | 2014 | 6 | 0.0812 | 0.0892 |
# vetor de datas
dates <- as.Date(paste(swc_top$Ano, swc_top$DOY),
format = "%Y %j")
swc_top$date <- dates
head(swc_top)| Ano | DOY | swc_11cm | swc_21cm | date |
|---|---|---|---|---|
| 2009 | 183 | 0.0959 | 0.0482 | 2009-07-02 |
| 2009 | 184 | 0.0931 | 0.0522 | 2009-07-03 |
| 2009 | 185 | 0.0893 | 0.0552 | 2009-07-04 |
| 2009 | 186 | 0.0855 | 0.0560 | 2009-07-05 |
| 2009 | 187 | 0.0817 | 0.0557 | 2009-07-06 |
| 2009 | 188 | 0.0785 | 0.0551 | 2009-07-07 |
Gerando série de datas de referência para combinar com dados de umidade do solo. Esse procedimento gerá uma série combinada com todos intervalos de tempo do período.
# seq dates de referencia
dates_ref <- data.frame(date = seq(min(dates), max(dates), "days"))
# combinando dados para o periodo de dates_ref
swc_top_all <- join(dates_ref, swc_top, by = "date", type = "left")Os dados de umidade do solo estão em m3 H2O/m3 solo. Vamos converter para lâmina de água equivalente.
# convertendo de umidade para armazenamento
swc_top_all$swctop <- rowMeans(swc_top_all[, c("swc_11cm","swc_21cm")])
swc_top_all$swctop <- swc_top_all$swctop*0.21*1000
summary(swc_top_all)| date | Ano | DOY | swc_11cm | swc_21cm | swctop | |
|---|---|---|---|---|---|---|
| Min. :2009-07-02 | Min. :2009 | Min. : 1 | Min. :0.0452 | Min. :0.0452 | Min. : 9.49 | |
| 1st Qu.:2010-10-01 | 1st Qu.:2010 | 1st Qu.: 92 | 1st Qu.:0.0656 | 1st Qu.:0.0607 | 1st Qu.:13.51 | |
| Median :2011-12-31 | Median :2011 | Median :183 | Median :0.0944 | Median :0.0929 | Median :19.64 | |
| Mean :2011-12-31 | Mean :2011 | Mean :183 | Mean :0.0977 | Mean :0.0953 | Mean :20.26 | |
| 3rd Qu.:2013-03-31 | 3rd Qu.:2013 | 3rd Qu.:274 | 3rd Qu.:0.1215 | 3rd Qu.:0.1205 | 3rd Qu.:24.91 | |
| Max. :2014-06-30 | Max. :2014 | Max. :366 | Max. :0.2392 | Max. :0.2391 | Max. :50.22 |
range(swc_top_all$date)[1] "2009-07-02" "2014-06-30" sum(is.na(swc_top_all$swctop))[1] 0Por prevenção vamos criar uma cópia dos dados swc_top_all chamada X.
# novo nome para os dados de swc
X <- swc_top_all
range(X$date)[1] "2009-07-02" "2014-06-30"head(X)| date | Ano | DOY | swc_11cm | swc_21cm | swctop |
|---|---|---|---|---|---|
| 2009-07-02 | 2009 | 183 | 0.0959 | 0.0482 | 15.1305 |
| 2009-07-03 | 2009 | 184 | 0.0931 | 0.0522 | 15.2565 |
| 2009-07-04 | 2009 | 185 | 0.0893 | 0.0552 | 15.1725 |
| 2009-07-05 | 2009 | 186 | 0.0855 | 0.0560 | 14.8575 |
| 2009-07-06 | 2009 | 187 | 0.0817 | 0.0557 | 14.4270 |
| 2009-07-07 | 2009 | 188 | 0.0785 | 0.0551 | 14.0280 |
E uma cópia dos dados nee_d.
# dados temporários
tmp <- nee_d
tmp$date <- as.Date(paste(nee_d$Year, nee_d$DoY), format = "%Y %j")
range(tmp$date)[1] "2009-07-08" "2014-07-01"head(tmp)| Year | DoY | Hour | NEE | PAR | Tair | VPD | Ustar | lai | swc1m | swc2m | date |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2009 | 189 | 50 | 2.11495 | 0 | 17.370 | 4.22681 | 0.580 | 1.5226 | 46.7653 | 123.205 | 2009-07-08 |
| 2009 | 189 | 100 | 4.38557 | 0 | 16.943 | 3.80730 | 0.645 | 1.5226 | 46.7653 | 123.205 | 2009-07-08 |
| 2009 | 189 | 150 | 2.11495 | 0 | 16.870 | 3.81223 | 0.618 | 1.5226 | 46.7653 | 123.205 | 2009-07-08 |
| 2009 | 189 | 200 | 2.40042 | 0 | 16.700 | 3.67781 | 0.533 | 1.5226 | 46.7653 | 123.205 | 2009-07-08 |
| 2009 | 189 | 250 | 1.73946 | 0 | 16.653 | 3.71807 | 0.520 | 1.5226 | 46.7653 | 123.205 | 2009-07-08 |
| 2009 | 189 | 300 | 3.13241 | 0 | 16.537 | 3.64160 | 0.504 | 1.5226 | 46.7653 | 123.205 | 2009-07-08 |
Combinando dados de 30 min com os diários, para repetir o armazenamento de água no topo do solo (camada superficial de 0 a 21 cm) para os 30 minutos de cada dia. Então selecionamos o período de dados de acordo com o conjunto de dados nee_d_swc, para que a série resultante (nee_d_swc) seja restrita ao mesmo período.
nee_d_swc <- join(tmp,
subset(X, sel = c("date", "swctop")),
by = "date",
type = "left")
nee_d_swc <- selectByDate(nee_d_swc,
start = format(max(c(min(tmp$date), min(X$date))), "%d/%m/%Y"),
end = format(min(c(max(tmp$date), max(X$date))), "%d/%m/%Y"))
# removendo date
nee_d_swc$date <- NULL
head(nee_d_swc)| Year | DoY | Hour | NEE | PAR | Tair | VPD | Ustar | lai | swc1m | swc2m | swctop |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2009 | 189 | 50 | 2.11495 | 0 | 17.370 | 4.22681 | 0.580 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 100 | 4.38557 | 0 | 16.943 | 3.80730 | 0.645 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 150 | 2.11495 | 0 | 16.870 | 3.81223 | 0.618 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 200 | 2.40042 | 0 | 16.700 | 3.67781 | 0.533 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 250 | 1.73946 | 0 | 16.653 | 3.71807 | 0.520 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 300 | 3.13241 | 0 | 16.537 | 3.64160 | 0.504 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
tail(nee_d_swc)| Year | DoY | Hour | NEE | PAR | Tair | VPD | Ustar | lai | swc1m | swc2m | swctop | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 87307 | 2014 | 181 | 2150 | 0.257 | 0.013667 | 17.5233 | 9.16804 | 0.107 | 1.0266 | 41.5677 | 88.6404 | 10.6575 |
| 87308 | 2014 | 181 | 2200 | -0.252 | 0.044667 | 17.2033 | 8.70899 | 0.107 | 1.0266 | 41.5677 | 88.6404 | 10.6575 |
| 87309 | 2014 | 181 | 2250 | -0.181 | 0.003333 | 17.0933 | 8.46399 | 0.042 | 1.0266 | 41.5677 | 88.6404 | 10.6575 |
| 87310 | 2014 | 181 | 2300 | 0.401 | 0.000000 | 17.0633 | 8.33440 | 0.053 | 1.0266 | 41.5677 | 88.6404 | 10.6575 |
| 87311 | 2014 | 181 | 2350 | NA | 0.000000 | 17.0000 | 8.35207 | 0.047 | 1.0266 | 41.5677 | 88.6404 | 10.6575 |
| 87312 | 2014 | 181 | 2400 | -0.054 | 0.000000 | 16.8067 | 8.12456 | 0.044 | 1.0266 | 41.5677 | 88.6404 | 10.6575 |
sum(is.na(nee_d_swc$swctop))/48[1] 05.4 Adição de coluna DateTime
Convertendo string de horas de “50, 100, 150,…” para “0.5, 1.0, 1.5, …” e adicionando time stamp no formato de tempo POSIX.
nee_d_swc$Hour <- change_hour_string(x = nee_d_swc$Hour)
head(nee_d_swc)| Year | DoY | Hour | NEE | PAR | Tair | VPD | Ustar | lai | swc1m | swc2m | swctop |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2009 | 189 | 0.5 | 2.11495 | 0 | 17.370 | 4.22681 | 0.580 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 1.0 | 4.38557 | 0 | 16.943 | 3.80730 | 0.645 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 1.5 | 2.11495 | 0 | 16.870 | 3.81223 | 0.618 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 2.0 | 2.40042 | 0 | 16.700 | 3.67781 | 0.533 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 2.5 | 1.73946 | 0 | 16.653 | 3.71807 | 0.520 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
| 2009 | 189 | 3.0 | 3.13241 | 0 | 16.537 | 3.64160 | 0.504 | 1.5226 | 46.7653 | 123.205 | 13.6185 |
nee_d_swc_px <- fConvertTimeToPosix(Data.F = nee_d_swc,
TFormat.s = 'YDH',
Year.s='Year',
Day.s='DoY',
Hour.s='Hour')5.5 Adicionando coluna com Radiaçao solar Global incidente
A radiação solar global incidente sera adicionada ao conjunto de dados usando a relação com PAR, obtida previamente como: \(R_g = 0.48PAR + 9.4 (PAR > 0, R^2=0.99)\).
nee_d_swc_px <- mutate(nee_d_swc_px,
Rg = ifelse(PAR > 0, 0.48*PAR + 9.4, PAR))
psych::describe(nee_d_swc_px[, c("PAR", "Rg")])| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| PAR | 1 | 87312 | 424.866 | 603.692 | 9.16733 | 312.538 | 13.6014 | -0.109 | 2680.33 | 2680.44 | 1.24408 | 0.298367 | 2.043046 |
| Rg | 2 | 87312 | 209.480 | 292.509 | 13.80032 | 155.775 | 20.4604 | -0.109 | 1295.96 | 1296.07 | 1.22890 | 0.259514 | 0.989925 |
5.6 Verificação da sequência cronológica dos dados
Verificação da série de dados de 30 min. A função fCheckHHTimeSeries() produz erros ou avisos caso algum problema de continuidade temporal na série seja detectado. Se estiver tudo ok ela não retorna nada.
fCheckHHTimeSeries(Time.V.p = nee_d_swc_px$DateTime,
DTS.n = 48,
CallFunction.s = "")5.7 Conversão dos dados para classe de dados do REddyProc
Para usar as funções do pacote REddyProc o conjunto de dados médios de 30 min devem ser convertidos para classe de referência (R5) denominada sEddyProc. O objeto criado (EddyProc.C) será usado como entrada nas funções disponíveis no pacote.
Ao inicializar uma classe sEddyProc, há uma séries de testes verificação sobre os dados, dos quais os mais relevantes são:
- Verificação das propriedades da série temporal (
?fCheckHHTimeSeries)
- Plausibilidade dos valores das variáveis (
?fCheckColPlausibility)
Vamos então gerar um objeto da classe sEddyProc a partir dos dados nee_d_swc_px e verificar se ocorre algum aviso relacionado as propriedades das séries ou a plausibilidade dos valores das variáveis.
EddyProc.C <- sEddyProc$new(ID.s = "PdG",
Data.F = nee_d_swc_px,
ColNames.V.s = c('NEE','PAR','Tair','VPD', 'Ustar', 'Rg'),
ColPOSIXTime.s = "DateTime")Warning in fCheckOutsideRange(Data.F,
VarName.V.s[v.i], c("<", 0), SubCallFunc.s):
sEddyProc.initialize:::fCheckColPlausibility:::fCheckOutsideRange:::
Variable outside (plausible) range in
1650 cases! Invalid values with 'PAR < 0':
-0.010,-0.011,-0.015,-0.005,-0.010,-0.005,-0.005,-0.005,-0.031,-0.046,-0.051,-0.025,-0.026,-0.061,-0.026,-0.005,-0.005,-0.005,-0.005,-0.005,-0.005,-0.005,-0.005,-0.020,-0.031,-0.005,-0.005,-0.005,-0.020,-0.046,-0.031,-0.056,-0.087,-0.062,-0.005,-0.046,-0.010,-0.010,-0.005,-0.005,-0.010,-0.005,-0.010,-0.010,-0.031,-0.005,-0.005,-0.003,-0.003,-0.010 ...Warning in fCheckOutsideRange(Data.F,
VarName.V.s[v.i], c(">", 2500), SubCallFunc.s):
sEddyProc.initialize:::fCheckColPlausibility:::fCheckOutsideRange:::
Variable outside (plausible) range in 5 cases! Invalid values with 'PAR >
2500': 2535,2564,2680,2542,2500 ...Warning in fCheckOutsideRange(Data.F,
VarName.V.s[v.i], c(">", 50), SubCallFunc.s):
sEddyProc.initialize:::fCheckColPlausibility:::fCheckOutsideRange:::
Variable outside (plausible) range in 2 cases! Invalid values with 'VPD >
50': 91,56 ...Warning in fCheckOutsideRange(Data.F,
VarName.V.s[v.i], c("<", 0), SubCallFunc.s):
sEddyProc.initialize:::fCheckColPlausibility:::fCheckOutsideRange:::
Variable outside (plausible) range in
1650 cases! Invalid values with 'Rg < 0':
-0.010,-0.011,-0.015,-0.005,-0.010,-0.005,-0.005,-0.005,-0.031,-0.046,-0.051,-0.025,-0.026,-0.061,-0.026,-0.005,-0.005,-0.005,-0.005,-0.005,-0.005,-0.005,-0.005,-0.020,-0.031,-0.005,-0.005,-0.005,-0.020,-0.046,-0.031,-0.056,-0.087,-0.062,-0.005,-0.046,-0.010,-0.010,-0.005,-0.005,-0.010,-0.005,-0.010,-0.010,-0.031,-0.005,-0.005,-0.003,-0.003,-0.010 ...Warning in fCheckOutsideRange(Data.F,
VarName.V.s[v.i], c(">", 1200), SubCallFunc.s):
sEddyProc.initialize:::fCheckColPlausibility:::fCheckOutsideRange:::
Variable outside (plausible) range in 7 cases! Invalid values with 'Rg >
1200': 1226,1240,1296,1230,1206,1207,1210 ...New sEddyProc class for site 'PdG'# estrutura da classe de dados
str(EddyProc.C$sDATA)'data.frame': 87312 obs. of 7 variables:
$ sDateTime: POSIXct, format: "2009-07-08 00:15:00" "2009-07-08 00:45:00" ...
$ NEE : num 2.11 4.39 2.11 2.4 1.74 ...
$ PAR : num 0 0 0 0 0 0 0 0 0 0 ...
$ Tair : num 17.4 16.9 16.9 16.7 16.7 ...
$ VPD : num 4.23 3.81 3.81 3.68 3.72 ...
$ Ustar : num 0.58 0.645 0.618 0.533 0.52 0.504 0.498 0.512 0.591 0.601 ...
$ Rg : num 0 0 0 0 0 0 0 0 0 0 ...Verificando algumas propriedades estatísticas das variáveis.
# resumo estatístico dos dados
round(describe(EddyProc.C$sDATA[,-1]), 2)| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NEE | 1 | 60173 | -2.15 | 8.03 | 0.20 | -1.57 | 5.28 | -49.88 | 39.80 | 89.69 | -0.56 | 2.17 | 0.03 |
| PAR | 2 | 87312 | 424.87 | 603.69 | 9.17 | 312.54 | 13.60 | -0.11 | 2680.33 | 2680.44 | 1.24 | 0.30 | 2.04 |
| Tair | 3 | 87312 | 22.29 | 4.58 | 21.96 | 22.29 | 4.59 | 4.66 | 36.62 | 31.96 | -0.02 | -0.16 | 0.02 |
| VPD | 4 | 87312 | 9.84 | 8.06 | 7.64 | 8.77 | 7.54 | 0.00 | 90.85 | 90.85 | 1.14 | 1.08 | 0.03 |
| Ustar | 5 | 76028 | 0.40 | 0.23 | 0.38 | 0.38 | 0.23 | 0.01 | 8.31 | 8.30 | 4.38 | 98.11 | 0.00 |
| Rg | 6 | 87312 | 209.48 | 292.51 | 13.80 | 155.77 | 20.46 | -0.11 | 1295.96 | 1296.07 | 1.23 | 0.26 | 0.99 |
5.8 Remoção ou Substituição de dados inconsistentes
Diante do aviso, vamos alterar os valores de PAR para as condições indicadas: PAR < 0 (1650) e PAR > 2500 (5 casos); VPD > 50 (2 casos); Rg > 1200 (7 casos); e reinicializar o processamento para o sítio do PdG.
# filtragem de PAR
nee_d_swc_px <- mutate(nee_d_swc_px,
PAR = pmin(PAR, 2500),
PAR = pmax(PAR, 0),
VPD = ifelse(VPD > 50, NA, VPD),
Rg = pmin(Rg, 1200),
Rg = pmax(Rg, 0))
# reinicialização
EddyProc.C <- sEddyProc$new(ID.s = "PdG",
Data.F = nee_d_swc_px,
ColNames.V.s = c('NEE','PAR','Rg','Tair','VPD', 'Ustar'),
ColPOSIXTime.s = "DateTime")New sEddyProc class for site 'PdG'# resumo estatístico dos dados
round(psych::describe(EddyProc.C$sDATA[, -1]), 2)| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NEE | 1 | 60173 | -2.15 | 8.03 | 0.20 | -1.57 | 5.28 | -49.88 | 39.80 | 89.69 | -0.56 | 2.17 | 0.03 |
| PAR | 2 | 87312 | 424.86 | 603.68 | 9.17 | 312.54 | 13.59 | 0.00 | 2500.00 | 2500.00 | 1.24 | 0.30 | 2.04 |
| Rg | 3 | 87312 | 209.48 | 292.50 | 13.80 | 155.77 | 20.46 | 0.00 | 1200.00 | 1200.00 | 1.23 | 0.26 | 0.99 |
| Tair | 4 | 87312 | 22.29 | 4.58 | 21.96 | 22.29 | 4.59 | 4.66 | 36.62 | 31.96 | -0.02 | -0.16 | 0.02 |
| VPD | 5 | 87310 | 9.84 | 8.05 | 7.64 | 8.77 | 7.54 | 0.00 | 47.31 | 47.31 | 1.12 | 0.97 | 0.03 |
| Ustar | 6 | 76028 | 0.40 | 0.23 | 0.38 | 0.38 | 0.23 | 0.01 | 8.31 | 8.30 | 4.38 | 98.11 | 0.00 |
Vamos reordenar as colunas de forma mais conveniente, para depois salvar os dados.
# reordenando colunas para que Rg e PAR fiquem próximas
nee_d_swc_px <- select(nee_d_swc_px, DateTime:PAR, Rg, Tair:lai, swctop, swc1m, swc2m)
head(nee_d_swc_px)| DateTime | Year | DoY | Hour | NEE | PAR | Rg | Tair | VPD | Ustar | lai | swctop | swc1m | swc2m |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2009-07-08 00:30:00 | 2009 | 189 | 0.5 | 2.11495 | 0 | 0 | 17.370 | 4.22681 | 0.580 | 1.5226 | 13.6185 | 46.7653 | 123.205 |
| 2009-07-08 01:00:00 | 2009 | 189 | 1.0 | 4.38557 | 0 | 0 | 16.943 | 3.80730 | 0.645 | 1.5226 | 13.6185 | 46.7653 | 123.205 |
| 2009-07-08 01:30:00 | 2009 | 189 | 1.5 | 2.11495 | 0 | 0 | 16.870 | 3.81223 | 0.618 | 1.5226 | 13.6185 | 46.7653 | 123.205 |
| 2009-07-08 02:00:00 | 2009 | 189 | 2.0 | 2.40042 | 0 | 0 | 16.700 | 3.67781 | 0.533 | 1.5226 | 13.6185 | 46.7653 | 123.205 |
| 2009-07-08 02:30:00 | 2009 | 189 | 2.5 | 1.73946 | 0 | 0 | 16.653 | 3.71807 | 0.520 | 1.5226 | 13.6185 | 46.7653 | 123.205 |
| 2009-07-08 03:00:00 | 2009 | 189 | 3.0 | 3.13241 | 0 | 0 | 16.537 | 3.64160 | 0.504 | 1.5226 | 13.6185 | 46.7653 | 123.205 |
6 Salvando dados de entrada para REddyProc
O conjunto de dados com médias de 30 min de NEE e variáveis ambientais dados nee_d_swc_px será salvo no formato de arquivo binário do R, cuja extensão do arquivo é .rds. Este formato permite carregar os dados facilmente no R.
saveRDS(object = nee_d_swc_px, file = "output/nee_pdg_reddypro_input.rds")7 Variação temporal das variáveis
Para usar as funções gráficas do pacote openair os dados devem ter uma coluna denominada date no formato POSIXt. Então vamos renomear a coluna DateTime do conjunto de dados nee_d_swc_px para date e fazer um gráfico das variáveis.
nee_d_swc_px_oa <- nee_d_swc_px
nee_d_swc_px_oa <- select(nee_d_swc_px_oa, -(Year:Hour))
nee_d_swc_px_oa <- rename(nee_d_swc_px_oa, date = DateTime)
timePlot(nee_d_swc_px_oa,
names(nee_d_swc_px_oa)[-1],
ylab = "variables",
date.format = "%b\n%Y",
date.breaks = 16,
date.pad = TRUE,
plot.type = "h",
#key.columns = 4,
key = FALSE)
Podemos visualizar as séries das variáveis para janelas de tamanho fixo, por exemplo N = 10 dias. Isso permite detectar padrões de variação das variáveis e identificar relações entre as variáveis.
N <- 48*10
interval <- (1:nrow(nee_d_swc_px_oa))%/%N
s <- split(nee_d_swc_px_oa, f = interval)
l <- l_ply(s,
function(x){
# x <- s[[1]]
tit <- paste(format(range(x$date), "%d %b %Y"),
collapse = " a ")
# gráficos dos N dias
timePlot(x,
names(x)[-1],
ylab = "variables",
date.format = "%d %b\n%Y",
#date.breaks = 16,
date.pad = TRUE,
key.columns = 4,
plot.type = "h",
main = tit)
# garbage clean
gc()
return(tit)
})