1 Introdução

Nesta etapa descreve-se a identificação de dados espúrios nas médias de 30 min. Foram criados para teste dois métodos:

Baseado na comparação entre a amplitude do ciclo diário mediano móvel (dentro de uma janela de N dias antes e após a meia-hora) e a diferença do valor da variável no instante atual e o seu valor no instante anterior. Se a diferença instantânea supera o valor da amplitude do ciclo diário uma flag de controle de qualidade é definida com o valor 1, se não seu valor é zero.
Baseado nos limiares superior e inferior do quantil da variável (p.ex. 1 e 99%) na meia-hora do dia. Se valor instantâneo da variável é menor (maior) que o quantil inferior (superior) uma flag de controle de qualidade é definida com valor 1, se não seu valor é zero.

Foram definidos dois cenários de filtragem para o NEE:

filtragem rigorosa: baseada na combinação do método do ciclo diário mediano móvel (com fator de amortecimento de 0.7 e a posterior aplicação do método dos quantis.
filtragem flexível: baseada no método dos quantis com limiar definido como 0.1 e 99.9 %.

Os dados espúrios não são eliminados da série, apenas são criados flags que indicam se o dados estão entre os limiares (0, dados de qualidade) ou não (1, dado espúrio).

Gráficos das séries das variáveis para cada cenário de filtragem foram gerados para facilitar a visualização do efeito das filtragens.

Por fim determinou-se a frequência de ocorrência de falhas e dados rejeitados nos dois cenários de filtragem e a distribuição das frequência entre os períodos diurnos e noturnos.

2 Pré-configuração: pacotes e scripts necessários

Carregando os pacotes necessários.

# data manipulation
pcks <- c("plyr", "dplyr","reshape2","lubridate","stringr",
          "openair","psych","REddyProc","ggplot2","scales",
          "lattice", "latticeExtra", "tidyr", "magrittr", 
          "dygraphs", "xts", "scales")
invisible(sapply(pcks, require, character.only = TRUE, quietly = TRUE))

Carregando scripts.

source("R/gaps.R") # find gaps length and start
source("R/to_oa.R") # convert to a open air object
source("R/regularize.R") # regularizes dates
source("R/qc_rcr.R")
source("R/qc_hh_quantile.R")
source("R/total_NA.R")
source("R/percent_NA.R")
source("R/auxs_funs.R")
source("R/plot_selected_dates.R")

3 Carregando dados

Vamos carregar os dados salvos na etapa anterior (1_prepDataNightNeeHourly.html) já num formato conveniente para conversão na classe de dados do pacote REddyProc (Reichstein and Moffat, 2014).

# loading night nee flux data pre-processed
head(nee_hh <- readRDS(file = "output/nee_pdg_reddypro_input.rds"))

DateTime	Year	DoY	Hour	NEE	Tair	VPD	Ustar	lai	swctop	swc1m	swc2m
2009-07-08 00:30:00	2009	189	0.5	2.114952	17.370	4.226812	0.580	1.5226	13.6185	46.76525	123.2052
2009-07-08 01:00:00	2009	189	1.0	4.385571	16.943	3.807304	0.645	1.5226	13.6185	46.76525	123.2052
2009-07-08 01:30:00	2009	189	1.5	2.114952	16.870	3.812228	0.618	1.5226	13.6185	46.76525	123.2052
2009-07-08 02:00:00	2009	189	2.0	2.400422	16.700	3.677811	0.533	1.5226	13.6185	46.76525	123.2052
2009-07-08 02:30:00	2009	189	2.5	1.739456	16.653	3.718069	0.520	1.5226	13.6185	46.76525	123.2052
2009-07-08 03:00:00	2009	189	3.0	3.132409	16.537	3.641602	0.504	1.5226	13.6185	46.76525	123.2052

# data inicial e final
range(nee_hh$DateTime)

[1] "2009-07-08 00:30:00 GMT" "2014-07-01 00:00:00 GMT"

Vamos visualizar a série de dados de 30 min do NEE.

# conversion to openair format
head(to_oa(nee_hh))

date	NEE	Tair	VPD	Ustar	lai	swctop	swc1m	swc2m
2009-07-08 00:30:00	2.114952	17.370	4.226812	0.580	1.5226	13.6185	46.76525	123.2052
2009-07-08 01:00:00	4.385571	16.943	3.807304	0.645	1.5226	13.6185	46.76525	123.2052
2009-07-08 01:30:00	2.114952	16.870	3.812228	0.618	1.5226	13.6185	46.76525	123.2052
2009-07-08 02:00:00	2.400422	16.700	3.677811	0.533	1.5226	13.6185	46.76525	123.2052
2009-07-08 02:30:00	1.739456	16.653	3.718069	0.520	1.5226	13.6185	46.76525	123.2052
2009-07-08 03:00:00	3.132409	16.537	3.641602	0.504	1.5226	13.6185	46.76525	123.2052

# regularized data
nee_hh_r <- regularize(to_oa(nee_hh))

Dataset time step: 30 mins

# gráfico da série temporal do NEE
timePlot(nee_hh_r, 
         "NEE",
         date.format = "%b\n%Y",
         date.breaks = 16, 
         date.pad = TRUE,         
         plot.type = "h")

4 Definiçao de períodos diurnos e noturnos

A radiação solar potencial incidente (PotRad) foi determinada pela função fCalcPotRadiation() do pacote REddyProc:

# coordenadas do PdG
lon <- -(47 + 38/60 + 2/3600)
lat <- -(21 + 37/60 + 16.1/3600)
# day of year
DoY.V.n <- as.numeric(format(nee_hh_r$date, '%j'))
# hour of day
Hour.V.n <- as.numeric(format(nee_hh_r$date, '%H')) + 
            as.numeric(format(nee_hh_r$date, '%M'))/60
# Potential Radiation
nee_hh_r$PotRad <- c(t(fCalcPotRadiation(DoY.V.n, 
                                     Hour.V.n, 
                                     Lat_deg.n = lat,
                                     Long_deg.n = lon,
                                     TimeZone_h.n = -3,
                                     useSolartime.b = TRUE)))

Gráfico de comparação entre a PotRad e Rg para Julho/2010 e Dezembro/2011.

# verificando concordancia da PotRad com PAR incidente
timePlot(selectByDate(nee_hh_r, 
                      month = 7, 
                      year = 2010,
                      day = 15:25), 
         c("Rg", "PotRad"), 
         group = TRUE, 
         lwd = 2, 
         lty =1)

timePlot(selectByDate(nee_hh_r, 
                      month = 12, 
                      year = 2011,
                      day = 15:25), 
         c("Rg", "PotRad"), 
         group = TRUE, 
         lwd = 2, 
         lty =1)

Gráfico de comparação entre PotRad e Rg para todo período.

tp_rad <- timePlot(nee_hh_r, 
                   c("Rg", "PotRad"), 
                   date.breaks = 16,
                   group = TRUE, 
                   date.pad = TRUE,
                   lwd = c(3,0.5),
                   lty = c(1,3),
                   date.format = "%b\n%Y")

A variável daytime será adicionada como uma nova coluna no data frame de dados para identificar se a meia-hora é dia ou noite.

O critério utilizado para definição de período diurno foi \(PotRad > 0\) determinada conforme Reichstein, Falge, Baldocchi, Papale, et al. (2005).

limiar_rad <- 0
nee_hh_r <- mutate(nee_hh_r,
             daytime = ifelse(PotRad > 0, #& Rg > limiar_rad,
             #daytime = ifelse(PotRad > 0 | Rg > limiar_rad,
                              "day", 
                              "night")
                   )
head(select(nee_hh_r, Rg, PotRad, daytime), 48)

Rg	PotRad	daytime
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
9.41248	0.00000	night
15.32080	33.74165	day
64.72000	176.89004	day
148.87216	313.80634	day
243.72016	442.14788	day
325.15984	559.71870	day
430.47184	664.50712	day
518.83984	754.72020	day
582.67984	828.81434	day
642.84016	885.52180	day
690.36016	923.87228	day
613.87984	943.20960	day
745.72000	943.20289	day
743.80000	923.85227	day
692.11984	885.48883	day
699.00016	828.76898	day
406.95184	754.66321	day
296.07184	664.43949	day
362.79184	559.64158	day
403.72000	442.06260	day
149.11216	313.71435	day
120.40816	176.79291	day
49.30720	33.64105	day
14.56480	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night
0.00000	0.00000	night

# timePlot(nee_hh_r %>% filter(daytime == "night"), "NEE", type = "hour")
# timePlot(nee_hh_r %>% filter(daytime == "day"), "NEE", type = "hour")
# nee_hh_r %>% filter(daytime == "day") %>% timeVariation("NEE")
# nee_hh_r %>% filter(daytime == "night") %>% timeVariation("NEE")

5 Controle de qualidade

Uma forma alternativa para visualizar dados espúrios (ou discrepantes) é o gráfico das séries temporais de cada meia-hora da variável (no exemplo abaixo o NEE).

timePlot(nee_hh_r, 
         "NEE", 
         type = "hour",
         date.format = "%b\n%Y",
         #date.breaks = 16, 
         date.pad = TRUE,         
         plot.type = "h",
         #y.relation = "free"
         layout = c(3,8)
         )

Foram criados dois critérios para indicar os dados espúrios que podem prejudicar o ajuste de modelos de regressão não linear (e.g. Falge, Baldocchi, Olson, Anthoni, et al. (2001) entre a variável dependente, NEE, e as variáveis independentes (PAR, Tair, VPD, etc).

5.1 Critérios para indicação de dados espúrios

O primeiro, determinado pela função qc_rcr(), baseia-se na comparação entre a diferença \(\delta NEE_{i}\) do NEE na meia-hora atual e a meia-hora anterior \(NEE_{i-1}\) com a amplitude do ciclo diário mediano móvel (\(\Delta NEE_{i}\)) com janela de N dias (p.ex.: 60 dias) em torno da meia-hora considerada. Se \(\delta NEE_{i} > \Delta NEE_{i}\) o valor da meia-hora é marcado como potencial dados espúrio. Mais detalhes do procedimento na seção 9.

O segundo, determinado pela função qc_hh_quantile(), baseia-se no quantil de NEE das meia-horas. Se o valor instantâneo é menor ou maior ao valor correspondente ao quantil de um limiar (p.ex. 1 e 99%, respectivamente) esse valor é indicado por um flag (1: dado espúrio, 0: dado válido). Mais detalhes do procedimento na seção 9.

Ambos critérios podem ser aplicados mais de uma vez.

5.1.1 Aplicação ao `NEE`

5.1.1.1 Filtro baseado na amplitude do ciclo diário

A função qc_rcr() foi aplicada aos dados de 30 min de NEE com uma janela móvel de 60 dias para cálculo do ciclo diário mediano. Esse valor alto foi escolhido baseado no tamanho da maior falha de dados, garantindo que existam dados para determinação do ciclo diário mediano.

# nee_hh_r filtrado uma vez
nee_hh_r_f <- qc_rcr(hhdata = select(nee_hh_r, date, NEE), 
                     window.size = 60,
                     mean.function = median,
                     sd.function = mad, 
                     weight.sd.function = 0,
                     aux.vars = TRUE,
                     damping = 1)
# nee_hh_r filtrado pela amplitude do ciclo usando damping = 0.7
nee_hh_r_f2 <- qc_rcr(hhdata = select(nee_hh_r, date, NEE),
                        #select(mutate(nee_hh_r_f, 
                        #                     NEE = ifelse(qc_rcr_NEE == 1, NA, NEE)), 
                        #              date, NEE), 
                     window.size = 60,
                     mean.function = median,
                     sd.function = mad, 
                     weight.sd.function = 0,
                     aux.vars = TRUE,
                     damping = 0.7)

# temp <- select(nee_hh_r_f, date, NEE, qc_rcr_NEE)
# temp$qc_rcr_NEE2 <- nee_hh_r_f2$qc_rcr_NEE
# timePlot(filter(temp, qc_rcr_NEE == 0), "NEE")
# timePlot(filter(temp, qc_rcr_NEE2 == 0), "NEE")

A saída da função qc_rcr(aux.vars = TRUE, damping = 1) é mostrada abaixo (apenas as primeiras linhas). Note que os dados são mantidos intactos e uma nova coluna com a variável qc_rcr_NEE é criada com valor binário (0 ou 1).

Nos gráficos mostrados os valores de NEE acima ou abaixo da banda sombreada são marcados como espúrios. Os flags marcados em vermelho identificam dados espúrios armazenados no data frame nee_hh_r_f2 gerado no trecho de comandos acima.

# dados de com flag gerada pelo filtro baseado na amplitude do ciclo diário
# e o parâmetro damping = 0.7 (0.7 * amplitude)
head(nee_hh_r_f2)

date	NEE	qc_upper	qc_lower
2009-07-08 00:30:00	2.114952	NA	NA
2009-07-08 01:00:00	4.385571	NA	NA
2009-07-08 01:30:00	2.114952	8.360264	-3.809898
2009-07-08 02:00:00	2.400422	8.360264	-4.504700
2009-07-08 02:30:00	1.739456	8.228840	-5.199502
2009-07-08 03:00:00	3.132409	8.097417	-4.943307

## dados para gráfico
data_plot <- nee_hh_r_f2 %>% 
  selectByDate(year = 2012) %>%
  mutate(qcf = ifelse(qc_rcr_NEE == 1, NEE, NA)) %>%
  select(-qc_rcr_NEE)
## dados no formato xts
data_plot <- as.xts(data_plot[, -1], 
                    order.by = data_plot[, 1])
# dynamic graph
dygraph(data_plot, main = "NEE médio 30 min") %>%
  dySeries(c("qc_lower", "NEE", "qc_upper")) %>% 
  dySeries("qcf", drawPoints = TRUE, pointSize = 4)%>%
  dyOptions(useDataTimezone = TRUE) %>%
  dyRangeSelector()%>% 
  dyAxis("y", label = "NEE") %>%
  dyOptions(colors = c("black","red"))

Note que o gráfico é interativo: movendo o mouse sobre a série os valores individuais são mostrados. É possível selecionar regiões do gráfico para dar zoom (duplo click para sair do zoom). Para uma janela de zoom selecionada, passando o mouse no centro da janela (aparece um cruz) é possível varrer a série com o tamanho da janela fixada.

5.1.1.2 Filtro baseado no quantil da meia-hora

O filtro baseado no quantil horário será aplicado nos dados de NEE brutos e para comparação com os dados de NEE marcados anteriormente pelo critério da amplitude do ciclo diário (usando damping = 0.7), ou seja, somente as linhas dos dados em que a flag qc_rcr_NEE = 0 nos dados nee_hh_r_f2. No gráfico do ciclo diário abaixo, valores acima (abaixo) da linha azul (vermelha) são marcados como espúrios.

# quantis p/ filtragem
q_inf_nee <- 0.001
q_sup_nee <- 0.999
# FILTRAGEM FLEXIVEL, pelo quantil para os dados brutos
nee_hh_r_ff <- qc_hh_quantile(X = nee_hh_r, 
                                      var.name = "NEE", 
                                      qtls = c(q_inf_nee, q_sup_nee), 
                                      plot.cycle = TRUE,
                                      suffix = "ff")

Dataset time step: 30 mins

head(nee_hh_r_ff)

date	NEE	qc_NEE_ff
2009-07-08 00:30:00	2.114952	0
2009-07-08 01:00:00	4.385571	0
2009-07-08 01:30:00	2.114952	0
2009-07-08 02:00:00	2.400422	0
2009-07-08 02:30:00	1.739456	0
2009-07-08 03:00:00	3.132409	0

# total de dados marcados
sum(nee_hh_r_ff$qc_NEE_ff == 1)

[1] 27355

# FILTRAGEM RIGOROSA (fr) 
# pelo quantil nos dados filtrados pelo mét. da amplit. 
# do ciclo com damping = 0.7
nee_hh_r_f2_reg <- regularize(filter(nee_hh_r_f2, qc_rcr_NEE == 0))

  time    freq
1   30 98.21 %
2   60  1.57 %
3   90  0.17 %
4  120  0.04 %
5  150  0.01 %
6  180     0 %
7  210     0 %
Dataset time step: 30 mins

# substituindo NAs por 1 na flag
nee_hh_r_f2_reg %<>% mutate(qc_rcr_NEE = ifelse(is.na(qc_rcr_NEE), 0, qc_rcr_NEE))
nee_hh_r_fr <- qc_hh_quantile(X = nee_hh_r_f2_reg, 
                                      var.name = "NEE", 
                                      qtls = c(q_inf_nee, q_sup_nee), 
                                      plot.cycle = TRUE,
                                      suffix = "fr")

Dataset time step: 30 mins

#rm(nee_hh_r_f2_reg)
head(nee_hh_r_fr)

date	NEE	qc_NEE_fr
2009-07-08 00:30:00	2.114952	0
2009-07-08 01:00:00	4.385571	0
2009-07-08 01:30:00	2.114952	0
2009-07-08 02:00:00	2.400422	0
2009-07-08 02:30:00	1.739456	0
2009-07-08 03:00:00	3.132409	0

# total de dados marcados
sum(nee_hh_r_fr$qc_NEE_fr == 1)

[1] 29125

A saída da função qc_hh_quantile() mostrada acima (apenas as primeiras linhas),os dados são mantidos intactos e uma nova coluna com a variável qc_NEE_suffix é criada com valores binários (0 ou 1).

Vamos criar um data frame unificado com as duas flags usadas para identificação de dados espúrios. As flags serão chamadas:

qc_nee_fr: controle de qualidade do NEE pela filtragem rigorosa, onde usou-se dois critérios de filtragem, primeiro o baseado na amplitude do ciclo diário mediano e sobre esses a filtragem pelo quantil da meia-hora;
qc_nee_ff: controle de qualidade do NEE pela filtragem flexível, onde usou-se apenas o critério do quantil da meia-hora;

#head(nee_hh_r_ff)
nee_hh_r_qc <- nee_hh_r %>% 
  bind_cols(qc_nee_fr = nee_hh_r_fr %>%
                        select(grep("qc", names(nee_hh_r_fr))),
            qc_nee_ff = nee_hh_r_ff %>%
                        select(grep("qc", names(nee_hh_r_ff)))) %>% 
  tbl_df
glimpse(nee_hh_r_qc)

Observations: 87312
Variables:
$ date      (time) 2009-07-08 00:30:00, 2009-07-08 01:00:00, 2009-07-0...
$ NEE       (dbl) 2.1149520, 4.3855706, 2.1149520, 2.4004219, 1.739455...
$ PAR       (dbl) 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.0...
$ Rg        (dbl) 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000...
$ Tair      (dbl) 17.370, 16.943, 16.870, 16.700, 16.653, 16.537, 16.4...
$ VPD       (dbl) 4.226812, 3.807304, 3.812228, 3.677811, 3.718069, 3....
$ Ustar     (dbl) 0.580, 0.645, 0.618, 0.533, 0.520, 0.504, 0.498, 0.5...
$ lai       (dbl) 1.5226, 1.5226, 1.5226, 1.5226, 1.5226, 1.5226, 1.52...
$ swctop    (dbl) 13.6185, 13.6185, 13.6185, 13.6185, 13.6185, 13.6185...
$ swc1m     (dbl) 46.76525, 46.76525, 46.76525, 46.76525, 46.76525, 46...
$ swc2m     (dbl) 123.2052, 123.2052, 123.2052, 123.2052, 123.2052, 12...
$ PotRad    (dbl) 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000...
$ daytime   (chr) "night", "night", "night", "night", "night", "night"...
$ qc_NEE_fr (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ qc_NEE_ff (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

5.1.2 Demais variáveis

Para as demais variáveis ambientais será aplicada a função qc_hh_quantile(). Para visualizar os potenciais dados espúrios são apresentados alguns gráficos de comparação usando dois limiares de quantis, referenciados também pelos sufixos:

filtragem flexível (ff): com percentis inferior e superior de 0.1 e 99.9, respectivamente.
filtragem rigorosa (fr): com percentis inferior e superior de 1 e 99, respectivamente.

Para aplicar a função qc_hh_quantile() a cada variável do data frame nee_hh_r_qc. Os armazenamentos de água no solo (swc*) e o índice de área foliar (lai) não foram obtidos na resolução temporal de 30 min. Eles são valores diários replicados para cada meia-hora, por isso a filtragem não se aplica.

Então identificação de dados espúrios será aplicada às variáveis ambientais: PAR, Rg, Tair, VPD e Ustar.

# definição de quantis p/ filtragem flexível
q_inf_vars_ff <- 0.001 
q_sup_vars_ff <- 0.999
# definição de quantis p/ filtragem rigorosa
q_inf_vars_fr <- 0.01
q_sup_vars_fr <- 0.99
# variáveis com dados de 30 min selecionadas
var_names <- select(nee_hh_r, PAR:Ustar) %>% names()
var_names

[1] "PAR"   "Rg"    "Tair"  "VPD"   "Ustar"

# looping nas variáveis para aplicar a função qc_hh_quantile()
l <- llply(var_names, 
           function(ivar) {
             cat(ivar, "\n")
              # ivar = "Rg"
              # dados de 30 min - filtragem rigorosa
               hhdata_fr <- qc_hh_quantile(X = nee_hh_r_qc, 
                                      var.name = ivar, 
                                      qtls = c(q_inf_vars_fr, q_sup_vars_fr), 
                                      plot.cycle = FALSE,
                                      suffix = "fr")[,-c(1,2)]
               # garbage collection
               gc()
               # dados de 30 min - filtragem flexível
               hhdata_ff <- qc_hh_quantile(X = nee_hh_r, 
                                      var.name = ivar, 
                                      qtls = c(q_inf_vars_ff, q_sup_vars_ff), 
                                      plot.cycle = FALSE,
                                      suffix = "ff")[,-c(1,2)]
               
               return(data.frame(hhdata_fr, hhdata_ff))
            }
      )

PAR 
Dataset time step: 30 mins 
Dataset time step: 30 mins 
Rg 
Dataset time step: 30 mins 
Dataset time step: 30 mins 
Tair 
Dataset time step: 30 mins 
Dataset time step: 30 mins 
VPD 
Dataset time step: 30 mins 
Dataset time step: 30 mins 
Ustar 
Dataset time step: 30 mins 
Dataset time step: 30 mins

As linhas indicam os quantis inferior (vermelho) e superior (azul) para os valores da variável em cada meia-hora do dia.

Dataset time step: 30 mins

Dataset time step: 30 mins

Dataset time step: 30 mins

Dataset time step: 30 mins

Dataset time step: 30 mins

Dataset time step: 30 mins

Dataset time step: 30 mins

Dataset time step: 30 mins

Dataset time step: 30 mins

Dataset time step: 30 mins

Precisamos combinar os dados da lista em um data frame com dados de todas as variáveis.

# combina nee_hh_r_qc e temp
nee_hh_r_qc <- tbl_df(data.frame(nee_hh_r_qc, bind_cols(l)))
# 1as linhas do dataframe com flags de qc 
glimpse(nee_hh_r_qc)

Observations: 87312
Variables:
$ date        (time) 2009-07-08 00:30:00, 2009-07-08 01:00:00, 2009-07...
$ NEE         (dbl) 2.1149520, 4.3855706, 2.1149520, 2.4004219, 1.7394...
$ PAR         (dbl) 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0...
$ Rg          (dbl) 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.000...
$ Tair        (dbl) 17.370, 16.943, 16.870, 16.700, 16.653, 16.537, 16...
$ VPD         (dbl) 4.226812, 3.807304, 3.812228, 3.677811, 3.718069, ...
$ Ustar       (dbl) 0.580, 0.645, 0.618, 0.533, 0.520, 0.504, 0.498, 0...
$ lai         (dbl) 1.5226, 1.5226, 1.5226, 1.5226, 1.5226, 1.5226, 1....
$ swctop      (dbl) 13.6185, 13.6185, 13.6185, 13.6185, 13.6185, 13.61...
$ swc1m       (dbl) 46.76525, 46.76525, 46.76525, 46.76525, 46.76525, ...
$ swc2m       (dbl) 123.2052, 123.2052, 123.2052, 123.2052, 123.2052, ...
$ PotRad      (dbl) 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.000...
$ daytime     (chr) "night", "night", "night", "night", "night", "nigh...
$ qc_NEE_fr   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_NEE_ff   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_PAR_fr   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_PAR_ff   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_Rg_fr    (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_Rg_ff    (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_Tair_fr  (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_Tair_ff  (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_VPD_fr   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_VPD_ff   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_Ustar_fr (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_Ustar_ff (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...

Vamos verificar a qualidade dos dados de radiação à noite, que pelo critério do quantil indicaram alguns casos de radiação acima de 600 \(W\: m^{-2}\) à noite. Para identificar dados inconsistentes à noite vamos alterar as flags de qc_Rg_fr e qc_PAR_fr para que indiquem dados espúrios quando PotRad = 0, hora >= 20 ou hora < 5 & Rg > 0.

# selecionando dados noturnos de Rg e PAR
s <- filter(nee_hh_r_qc, qc_Rg_fr == 0) %>%
     filter((hour(date) >= 20 | hour(date) < 5 ) & PotRad == 0 & PAR > 0) %>%
     select(date, Rg, PAR, PotRad)
slice(s, which.max(s$PAR))

date	Rg	PAR	PotRad
2010-06-26 20:00:00	9.47536	0.157	0

nee_hh_r_qc %>%
 selectByDate(start = "29/01/2010", 
              end = "05/02/2010") %>% 
  timePlot(c("PAR", "PotRad"), date.pad = TRUE)

filter(nee_hh_r_qc, qc_Rg_fr == 0) %>%
 selectByDate(start = "29/01/2010", 
              end = "05/02/2010") %>% 
  timePlot(c("PAR", "PotRad"), date.pad = TRUE)

# filtragem de Rg noturno incluída no qc_Rg_fr e qc_PAR_fr
nee_hh_r_qc %<>% 
  mutate(qc_PAR_fr = ifelse((hour(date) > 20 | 
                            hour(date) < 5 ) & 
                            (PotRad == 0 & PAR > 0),
                            1, qc_PAR_fr))
timePlot(filter(nee_hh_r_qc, qc_PAR_fr == 0), 
         "PAR", 
         type = "hour",
         date.format = "%m\n%y",
         y.relation = "free",
         layout = c(4,6)
         )

timePlot(filter(nee_hh_r_qc, qc_PAR_fr == 0), 
         "PAR", 
         type = "hour",
         date.format = "%m\n%y",
         y.relation = "same",
         layout = c(4,6))

Finalmente, salvamos o data frame combinado que inclui as flags para os diferentes critérios.

saveRDS(nee_hh_r_qc, file = "output/nee_hh_r_qc_pdg.rds")

5.1.3 Série temporal de `NEE` de 30 min com indicação de dados espúrios

Para visualizar o efeito de usarmos uma filtragem flexível (ff, quantis de 0.1 e 99.9 %) ou rigorosa (fr, após filtragem pela amplitude ciclo diário) na série de NEE são mostrados os gráficos das séries temporais para cada caso.

Análise visual dos dados NEE identificados como espúrios pela função qc_rcr().

# range NEE bruto
rng <- range(select(nee_hh_r_qc, NEE), na.rm = TRUE)
# gráfico do NEE com filtragem pela amplitude do ciclo diário
timePlot(filter(nee_hh_r_qc, qc_NEE_ff == 0), 
         "NEE",
         date.pad = TRUE,
         ylim = rng)

# gráfico do NEE com filtragem pela amplitude do ciclo diário e quantil de 0.1%
timePlot(filter(nee_hh_r_qc, qc_NEE_fr == 0), 
         "NEE",
         date.pad = TRUE,
         ylim = rng)

5.1.4 Série temporal de `NEE` por hora.

O efeito da filtragem pode ser mais evidente ao visualizarmos a o NEE médio de 30 min separadamente para cada hora do dia.

timePlot(filter(nee_hh_r_qc, qc_NEE_ff == 0), 
         "NEE", 
         type = "hour",
         date.format = "%b\n%Y",
         #date.breaks = 16, 
         date.pad = TRUE,         
         plot.type = "h",
         #y.relation = "free"
         ylim = rng,
         layout = c(3,8)
         )

# range dos dados nee_hh_r_filt1, para definição do ylim 
# do gráfico dos dados nee_hh_r_filt2
timePlot(filter(nee_hh_r_qc, qc_NEE_fr == 0), 
         "NEE", 
         type = "hour",
         date.format = "%m\n%Y",
         #date.breaks = 16, 
         date.pad = TRUE,         
         plot.type = "h",
         #y.relation = "free"
         ylim = rng,
         layout = c(3,8)
         )

5.1.5 Série temporal de `NEE` médio noturno para os critérios de filtragem

5.1.5.1 Gráfico da Frequência de NEE noturno disponível

Vamos determinar a frequência de ocorrência de dados noturnos de NEE, para avaliar o peso da média do NEE noturno.

Começamos com o data frame com dados noturnos.

# data frame de noites
nights <- nee_hh_r_qc %>% 
  select(date, NEE, daytime, matches("qc_NEE")) %>%
  mutate(date_ymd = as.Date(date)) %>%
  filter(daytime == "night")

Agrupamos o data frame por dia e calculamos as frequências de ocorrência de dados de NEE válidos (não faltantes), dados que atenderam a filtragem rigorosa, a flexível e as médias de NEE para as respectivas filtragens.

# data frame de noites agrupado por data
nights_dates <- 
nights %>% 
  group_by(date_ymd) %>%
  summarise(# num. de intervalos de 30 min noturnos
            n = n(),
            # frequência de dados válidos
            freq_valid = sum(!is.na(NEE), na.rm = TRUE),
            # frequência de dados válidos pela fr
            freq_qc_fr = sum(qc_NEE_fr == 0),
            # frequência de dados válidos pela ff
            freq_qc_ff = sum(qc_NEE_ff == 0),
            # NEE médio para os dados válidos pela fr
            NEE_avg_fr = mean(NEE[qc_NEE_fr == 0]),
            # NEE médio para os dados válidos pela ff
            NEE_avg_ff = mean(NEE[qc_NEE_ff == 0])) %>%
  mutate(date = as.POSIXct(paste(date_ymd, "23:00:00"), tz = "GMT"),
         date_ymd = NULL)
# cols de interesse
cols <- c("date", names(nights_dates)[-ncol(nights_dates)])
# ordenando colunas
nights_dates <- select(nights_dates, one_of(cols))

Gráficos das frequências noturnas de dados por dia.

ss_freq <- select(nights_dates, -matches("NEE_avg"))
timePlot(ss_freq, 
         names(ss_freq)[-1], 
         date.format = "%b\n%Y",
         date.breaks = 18, 
         date.pad = TRUE,
         #group = TRUE,
         lwd = c(2, 2),
         lty = 1,
         ref.y = list(h = 24, col = 5, lty = 3, lwd = 3)
         )

A frequência de dados dados válidos por noite fica menos intermitente a partir de 2012, apesar de conter falhas de dados extensas.

O nindica o nº de meia-horas por noite e depende do critério usado para definir período diurno e noturno.

Vamos verificar a frequência de ocorrência de n (n = 24 corresponde a 12 horas noturnas).

barplot(table(nights_dates$n))

5.1.5.2 Gráficos das médias noturnas de NEE por dia.

ss_avg <- select(nights_dates, date, matches("NEE"))
timePlot(ss_avg, 
         names(ss_avg)[-1], 
         date.format = "%b\n%Y",
         date.breaks = 18, 
         date.pad = TRUE,
         ylim = range(ss_avg$NEE_avg_ff, na.rm = TRUE),
         lwd = c(2, 2),
         lty = 1,
         ref.y = list(h = 0, col = 5, lty = 3, lwd = 2)
         )

Pelo critério fr ocorreram médias noturnas de NEE < 0. Vamos navegar por essas datas no gráfico da série de NEE.

# datas com NEE < 0
dates_nee_fr_lt0 <- nights_dates %>% 
                    filter(NEE_avg_fr < 0) %>% 
                    select(date)
(dates_nee_fr_lt0 <- as.Date(c(t(dates_nee_fr_lt0))))

 [1] "2009-12-31" "2010-01-17" "2010-03-31" "2010-09-15" "2011-01-26"
 [6] "2011-02-02" "2011-04-06" "2011-04-14" "2011-05-07" "2012-03-24"
[11] "2013-01-17" "2013-05-30" "2014-07-01"

## dados para dygraph
data_plot <- transmute(nee_hh_r_qc, 
                       date = date, 
                       NEE = ifelse(qc_NEE_fr == 1, NA, NEE))
data_plot <- as.xts(x = data_plot[, -1], 
                    order.by = data_plot[[1]])
# dynamic graph
dygraph(data_plot, main = "NEE médio 30 min") %>%
  dySeries("NEE") %>%   
  dyOptions(useDataTimezone = TRUE) %>% 
  dyRangeSelector() %>%
  dyAxis("y", label = "NEE") %>%
  dyLimit(limit = 0, color = "red")

Navegando entre 11-21 de julho e 2012 observa-se 10 noites consecutivas com valores de NEE idênticos. Por que? Maioria dos casos em que a média do NEE noturno foi negativo foi devido a pouca disponibilidade de dados na noite para composição da média

Pelo critério fr ocorreram alguns casos nas médias noturnas que NEE > 10. Vamos encontrar essas datas.

# datas com NEE > 10
dates_nee_fr_gt0 <- nights_dates %>% 
                    filter(NEE_avg_fr > 10) %>% 
                    select(date)
(dates_nee_fr_gt0 <- as.Date(c(t(dates_nee_fr_gt0))))

 [1] "2009-11-29" "2009-12-03" "2009-12-08" "2009-12-13" "2009-12-14"
 [6] "2010-01-14" "2011-01-11" "2011-01-12" "2011-02-26" "2011-12-15"
[11] "2012-01-01" "2012-01-14" "2012-01-15" "2012-12-11" "2013-04-05"
[16] "2013-12-10" "2013-12-11" "2013-12-28" "2014-02-25"

## plots
plot_selected_dates(hh_data = mutate(nee_hh_r_qc, 
                                     date = date, 
                                     NEE = ifelse(qc_NEE_fr == 1, NA, NEE)), 
                    dts = dates_nee_fr_gt0,
                    vars = c("NEE","PAR","Rg","PotRad"), 
                    window = 1)

Maioria dos casos em que a média do NEE noturno passou de 10 \(\mu mol\: CO_{2}\: m^{2} s^{-1}\) foi devido a pouca de quantidades de dados noturnos na composição da média

6 Frequência de ocorrência de dados faltantes e rejeitados

Cálculo da frequência de ocorrência (absoluta e relativa em %) de dados faltantes e rejeitados pelos critérios de filtragem.

Procedimento para os dados faltantes.

freqNA <- nee_hh_r_qc %>% 
  select(NEE, one_of(var_names)) %>%
  summarise_each(funs(freq.abs = total_NA, 
                      freq.rel = percent_NA), 
                 everything()) %>%
  gather(variable, n) %>%
  mutate(missing = round(n, 2), n = NULL)

Procedimento para os dados rejeitados.

# dados rejeitados
freq_table_outl <- nee_hh_r_qc %>%  
  select(contains("qc", ignore.case = TRUE)) %>%
  gather(variable, flag) %>% 
  mutate(variable = gsub("qc_", "qc", variable)) %>%
  group_by(variable) %>%
  summarise(absolute = sum(flag==1), 
            relative = (sum(flag==1)/n())*100) %>%
  separate(col = variable, into = c("variable", "quantil"), sep = "_") %>%
  gather(key = frequency, 
         value = value, 
         absolute:relative) %>% 
  mutate(value = round(value, 2)) %>%
  spread(key = quantil, 
         value = value) %>% 
  arrange(frequency)

Combinação das duas tabelas de frequência de ocorrência.

# tabela combinada
tab_freq <- bind_cols(freq_table_outl, freqNA[,"missing"]) %>%
                select(variable, frequency, missing, ff, fr)

Tabela 5. Frequência absoluta de dados rejeitados (e faltantes) para os critérios de filtragem flexível (quantil de 0.1) e rigorosa (amplitude diária e quantil).
variable	missing	ff	fr
qcNEE	27139	27355	29125
qcPAR	0	2508	5899
qcRg	0	2478	3230
qcTair	0	187	966
qcUstar	2	11451	11928
qcVPD	11284	209	967

Tabela 6. Frequência relativa (%) de dados rejeitados (e faltantes) para os critérios de filtragem flexível e rigorosa.
variable	missing	ff	fr
qcNEE	31.08	31.33	33.36
qcPAR	0.00	2.87	6.76
qcRg	0.00	2.84	3.70
qcTair	0.00	0.21	1.11
qcUstar	0.00	13.12	13.66
qcVPD	12.92	0.24	1.11

A porcentagem de dados faltantes na série de 30 min do NEE do sítio PdG é de 31 % (ou 27139 períodos de 30 min) de um total de 87312 períodos de 30 min (~1819 dias).

A frequência de dados rejeitados e faltantes no caso de uma filtragem mais rigorosa do NEE é de apenas 33.4 %, ou seja uma perda de dados de apenas 2.28% (1986 intervalos de 30 min) em relação a frequência de dados faltantes.

Vetor com nome das variáveis que apresentaram dados faltantes.

(var_miss_data <- unique(c(t(
  tab_freq %>% 
  filter(missing > 0) %>%
  select(variable) %>%
  mutate(variable = gsub("qc","",variable))  
))))

[1] "NEE"   "Ustar" "VPD"

6.1 Dia e Noite

Cálculo da Frequência absoluta e relativa (%) de dados faltantes e dados rejeitados pelos diferentes critérios de filtragem para os turnos diurno e noturno.

Entre as variáveis para calcular a frequência de dados rejeitados vamos incluir PAR que pode vir a ser usada para o preenchimento de falhas no período diurno. Os resultados para essa variável são válidos também para Rg, uma vez que Rg foi derivada de PAR.

# variaveis selecionadas
sel_vars <- c(var_miss_data, "PAR")

Procedimento para o cálculo de dados faltantes por turno.

# colunas de interesse
cols <- names(tbl_df(nee_hh_r))[-c(1, ncol(nee_hh_r))]
# dados faltantes
freqNA_dn <- tbl_df(nee_hh_r) %>%  
  gather(variable, value, one_of(cols)) %>% 
  mutate(variable = as.character(variable)) %>%
  filter(variable %in% sel_vars) %>%
  group_by(variable) %>%
  mutate(n_miss = total_NA(value), 
         n_data = n()) %>%
  group_by(variable, daytime) %>%
  summarise(freq.abs = (total_NA(value)), 
            n.miss.all = unique(n_miss), 
            n = unique(n_data)) %>%
  mutate(freq.rel = (freq.abs/n.miss.all)*100, 
         freq.rel.all = (freq.abs/n)*100) %>%
  select(-c(n, n.miss.all)) %>%
  gather(frequency, missing, freq.abs:freq.rel.all) %>%
  mutate(frequency = gsub("freq.abs","absolute",frequency),
         frequency = gsub("freq.rel","relative",frequency),
         missing = round(missing, 2),
         missing = ifelse(is.na(missing), 0, missing)) %>%
  arrange(frequency, variable)

Procedimento para o cálculo de dados rejeitados por turno.

# dados rejeitados
freq_table_outl_dn <- 
  nee_hh_r_qc %>%  
  select(daytime, contains("qc", ignore.case = TRUE)) %>%
  gather(variable, flag, contains("qc", ignore.case = TRUE)) %>% 
    mutate(variable = gsub("qc_", "qc", variable)) %>%
     group_by(variable) %>%
       mutate(n_rej = sum(flag==1), 
              n_data = n()) %>%
        group_by(variable, daytime) %>%
         summarise(freq.abs = sum(flag==1),
                   n.rej.all = unique(n_rej),
                   n = unique(n_data)) %>%
          separate(col = variable, 
                   into = c("variable", "quantil"), 
                   sep = "_") %>%
           mutate(freq.rel = (freq.abs/n.rej.all)*100, 
                  freq.rel.all = (freq.abs/n)*100) %>%
            select(-c(n, n.rej.all)) %>%
             gather(key = frequency, 
                    value = value, 
                    freq.abs:freq.rel.all) %>%
             mutate(value = round(value, 2)) %>%
              spread(key = quantil, 
                     value = value) %>% 
                mutate(frequency = gsub("freq.abs", "absolute", frequency),
                       frequency = gsub("freq.rel", "relative", frequency)) %>%
                 filter(variable %in% paste0("qc", sel_vars)) %>%
                 arrange(frequency, variable)

Combinação das duas tabelas de frequência de ocorrência.

# tabela combinada com colunas reordenadas
tab_freq_dn <- bind_cols(freq_table_outl_dn, freqNA_dn[,"missing"]) %>%
                select(variable, daytime, frequency, missing, ff, fr)

Tabela 7. Frequência absoluta de dados rejeitados (e faltantes) em relação ao total de dados rejeitados, para os critérios de filtragem flexível e rigorosa por turno (dia e noite).
variable	daytime	missing	ff	fr
qcNEE	day	12616	12721	13714
qcNEE	night	14523	14634	15411
qcPAR	day	0	2285	2728
qcPAR	night	0	223	3171
qcUstar	day	5712	5798	6069
qcUstar	night	5572	5653	5859
qcVPD	day	2	159	640
qcVPD	night	0	50	327

Tabela 8. Frequência relativa de dados rejeitados (e faltantes) em relação ao total de dados rejeitados, para os critérios de filtragem flexível e rigorosa por turno (dia e noite).
variable	daytime	missing	ff	fr
qcNEE	day	46.49	46.50	47.09
qcNEE	night	53.51	53.50	52.91
qcPAR	day	0.00	91.11	46.25
qcPAR	night	0.00	8.89	53.75
qcUstar	day	50.62	50.63	50.88
qcUstar	night	49.38	49.37	49.12
qcVPD	day	100.00	76.08	66.18
qcVPD	night	0.00	23.92	33.82

Tabela 9. Frequência de dados rejeitados (e faltantes) em relação ao total de observações da variável, para os critérios de filtragem flexível e rigorosa por turno (dia e noite).
variable	daytime	missing	ff	fr
qcNEE	day	14.45	14.57	15.71
qcNEE	night	16.63	16.76	17.65
qcPAR	day	0.00	2.62	3.12
qcPAR	night	0.00	0.26	3.63
qcUstar	day	6.54	6.64	6.95
qcUstar	night	6.38	6.47	6.71
qcVPD	day	0.00	0.18	0.73
qcVPD	night	0.00	0.06	0.37

Há um amento na frequência de dados rejeitados no período noturno (diurno) de ~ 1 % (1 %), em relação a frequência e dados faltantes pela filtragem rigorosa.

Com as informações de frequência de ocorrência de falhas e rejeições para os diferentes critérios; em combinação com os gráficos das variáveis é possível optar por um dos critérios buscando o equilíbrio entre o rigor na filtragem (menos dados espúrios) e o aumento de dados desprezados.

7 Tamanho das falhas

A posição e o tamanho das falhas de nee_obs são identificadas pela data inicial (sdate) e tamanho (length), respectivamente na tabela abaixo (suprimida).

## tamanho das falhas de NEE
g <- nee_hh_r_qc %$% gaps(NEE)
## adicionando data de início da falha
g$sdate <- nee_hh_r_qc %$% date[g$start]
## ordenando pelo tamanho da falha
gaps_tab <- g %$% length %>% 
            order(decreasing = TRUE) %>%
            slice(g, .)
## adicionando coluna de fim da falha
gaps_tab %<>% 
  mutate(sdate = as.Date(sdate),
         edate = sdate + length - 1)
kable(head(gaps_tab, 20), 
      row.names = TRUE, 
      caption = "Tabela 10. Data e tamanho das falhas de NEE.  A primeira coluna da tabela indica o tamanho da falha.")

Tabela 10. Data e tamanho das falhas de NEE. A primeira coluna da tabela indica o tamanho da falha.
	length	start	sdate	edate
1	2010	56482	2012-09-26	2018-03-28
2	1318	75615	2013-10-30	2017-06-08
3	1306	9814	2010-01-28	2013-08-25
4	1254	50424	2012-05-23	2015-10-28
5	683	20855	2010-09-15	2012-07-28
6	647	13701	2010-04-19	2012-01-25
7	596	18930	2010-08-06	2012-03-23
8	535	63051	2013-02-10	2014-07-29
9	486	7999	2009-12-21	2011-04-20
10	466	38774	2011-09-23	2012-12-31
11	437	15677	2010-05-30	2011-08-09
12	360	45142	2012-02-03	2013-01-27
13	237	6450	2009-11-19	2010-07-13
14	153	17447	2010-07-06	2010-12-05
15	153	82645	2014-03-25	2014-08-24
16	148	4974	2009-10-19	2010-03-15
17	143	13002	2010-04-04	2010-08-24
18	126	39865	2011-10-16	2012-02-18
19	117	40502	2011-10-29	2012-02-22
20	111	33364	2011-06-03	2011-09-21

Gráfico da porcentagem de falhas pelo tamanho das falhas.

# tabela de num. falhas por tamanho
gfreqlen <- g %$% 
  table(g$length) %>% 
  #prop.table() %>% `*`(100) %>%
  as.data.frame(stringsAsFactors = FALSE) %>%
  rename(gap_lenght = Var1, n_gaps = Freq) %>%
  mutate(gap_lenght = as.integer(as.character(gap_lenght))) 
# gráfico do tamanho pelo num. de falhas
gfreqlen %>% 
  ggplot(aes(gap_lenght, n_gaps)) + 
  geom_line() + geom_point(aes(size = 2)) +
  scale_y_continuous(breaks = pretty_breaks(10)) + 
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  theme_bw(base_size = 15) +
   theme(legend.position="none") + 
  labs(x = "Tamanho da falha (meia-horas)", y = "Nº de falhas")

# tabela de num. falhas (%) por tamanho
gfreqlen_pctg <- g %$% 
  table(g$length) %>% 
  prop.table() %>% `*`(100) %>%
  as.data.frame(stringsAsFactors = FALSE) %>%
  rename(gap_lenght = Var1, n_gaps = Freq) %>%
  mutate(gap_lenght = as.integer(as.character(gap_lenght))) 
# gráfico do tamanho (log10) pelo num. de falhas
gfreqlen_pctg %>% 
  ggplot(aes(gap_lenght, n_gaps)) + 
  geom_line() + geom_point(aes(size = 2)) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x)),
                limits = c(0.03, 100)) +
  theme_bw(base_size = 15) +
  theme(legend.position="none") + 
  labs(x = "Tamanho da falha (meia-horas)", y = "Nº de falhas (%)")

Determinação da frequência de ocorrência das classes de tamanho das falhas por ano.

freq_year <- g %$% 
             table(length, year(sdate)) %>% 
             addmargins(margin = 2)
kable(freq_year, 
      row.names = TRUE, 
      caption = "Tabela 11. Frequência de ocorrência de falhas por classes de tamanho das falhas e ano.  A primeira coluna da tabela indica o tamanho da falha.")

Tabela 11. Frequência de ocorrência de falhas por classes de tamanho das falhas e ano. A primeira coluna da tabela indica o tamanho da falha.
	2009	2010	2011	2012	2013	2014	Sum
1	121	303	417	131	98	38	1108
2	80	139	222	70	54	26	591
3	24	87	122	36	37	9	315
4	19	51	81	19	20	7	197
5	17	40	46	12	16	8	139
6	17	36	45	4	13	8	123
7	17	22	31	11	14	5	100
8	8	16	14	2	9	2	51
9	4	15	26	3	8	5	61
10	7	13	9	2	8	2	41
11	5	7	8	3	5	5	33
12	3	5	13	0	1	1	23
13	7	4	12	1	4	4	32
14	4	11	9	1	2	0	27
15	2	5	11	0	2	1	21
16	0	5	6	1	2	3	17
17	1	4	5	1	0	1	12
18	2	2	6	0	4	0	14
19	2	3	4	1	3	1	14
20	3	3	4	0	4	0	14
21	1	2	6	2	1	0	12
22	0	4	1	0	2	0	7
23	2	2	5	2	0	0	11
24	1	0	3	0	1	0	5
25	2	3	1	1	0	0	7
26	0	2	4	0	0	1	7
27	1	4	3	1	0	0	9
28	1	3	3	0	0	0	7
29	0	2	5	1	1	0	9
30	1	1	0	0	2	0	4
31	0	1	2	0	1	0	4
32	0	3	2	0	0	0	5
33	0	3	2	1	1	1	8
34	1	2	0	0	0	0	3
35	0	1	1	0	0	0	2
36	1	1	1	0	1	0	4
37	1	1	1	0	0	0	3
38	0	1	0	0	1	0	2
39	0	0	2	0	0	0	2
40	1	1	3	0	1	0	6
41	1	1	1	0	0	0	3
42	1	0	2	0	1	0	4
43	0	1	1	0	0	0	2
45	0	0	1	1	0	0	2
46	0	0	1	0	0	0	1
50	0	1	1	0	0	0	2
51	2	0	0	0	0	0	2
52	1	1	0	0	0	0	2
55	1	0	0	0	0	0	1
58	1	0	1	0	0	0	2
61	1	0	0	0	0	0	1
62	0	1	0	0	1	0	2
63	0	0	1	0	0	0	1
68	0	0	1	0	0	0	1
71	0	0	1	0	0	0	1
82	0	1	2	0	0	0	3
83	0	0	1	0	0	0	1
84	0	0	1	0	0	0	1
96	0	0	1	0	0	0	1
97	0	1	0	0	0	0	1
111	0	0	1	0	0	0	1
117	0	0	1	0	0	0	1
126	0	0	1	0	0	0	1
143	0	1	0	0	0	0	1
148	1	0	0	0	0	0	1
153	0	1	0	0	0	1	2
237	1	0	0	0	0	0	1
360	0	0	0	1	0	0	1
437	0	1	0	0	0	0	1
466	0	0	1	0	0	0	1
486	1	0	0	0	0	0	1
535	0	0	0	0	1	0	1
596	0	1	0	0	0	0	1
647	0	1	0	0	0	0	1
683	0	1	0	0	0	0	1
1254	0	0	0	1	0	0	1
1306	0	1	0	0	0	0	1
1318	0	0	0	0	1	0	1
2010	0	0	0	1	0	0	1

Determinação da frequência de ocorrência do total de falhas por ano.

# frequencia anual
freq_tot_year <- g %$% 
                 table(length, year(sdate)) %>%
                 as.matrix() 
# adicionando margins
freq_tot_year_margins <- freq_tot_year %>% 
  rownames() %>% 
   as.integer() * freq_tot_year %>%
   cbind(., freq_tot_year %>% rowSums()) 
freq_tot_year_margins %<>% colSums() %>% rbind(freq_tot_year_margins, .)
# tabela
kable(freq_tot_year, 
      row.names = TRUE, 
      caption = "Tabela 12. Total de falhas por classes de tamanho das falhas e ano.  A primeira coluna da tabela indica o tamanho da falha.")

Tabela 12. Total de falhas por classes de tamanho das falhas e ano. A primeira coluna da tabela indica o tamanho da falha.
	2009	2010	2011	2012	2013	2014
1	121	303	417	131	98	38
2	80	139	222	70	54	26
3	24	87	122	36	37	9
4	19	51	81	19	20	7
5	17	40	46	12	16	8
6	17	36	45	4	13	8
7	17	22	31	11	14	5
8	8	16	14	2	9	2
9	4	15	26	3	8	5
10	7	13	9	2	8	2
11	5	7	8	3	5	5
12	3	5	13	0	1	1
13	7	4	12	1	4	4
14	4	11	9	1	2	0
15	2	5	11	0	2	1
16	0	5	6	1	2	3
17	1	4	5	1	0	1
18	2	2	6	0	4	0
19	2	3	4	1	3	1
20	3	3	4	0	4	0
21	1	2	6	2	1	0
22	0	4	1	0	2	0
23	2	2	5	2	0	0
24	1	0	3	0	1	0
25	2	3	1	1	0	0
26	0	2	4	0	0	1
27	1	4	3	1	0	0
28	1	3	3	0	0	0
29	0	2	5	1	1	0
30	1	1	0	0	2	0
31	0	1	2	0	1	0
32	0	3	2	0	0	0
33	0	3	2	1	1	1
34	1	2	0	0	0	0
35	0	1	1	0	0	0
36	1	1	1	0	1	0
37	1	1	1	0	0	0
38	0	1	0	0	1	0
39	0	0	2	0	0	0
40	1	1	3	0	1	0
41	1	1	1	0	0	0
42	1	0	2	0	1	0
43	0	1	1	0	0	0
45	0	0	1	1	0	0
46	0	0	1	0	0	0
50	0	1	1	0	0	0
51	2	0	0	0	0	0
52	1	1	0	0	0	0
55	1	0	0	0	0	0
58	1	0	1	0	0	0
61	1	0	0	0	0	0
62	0	1	0	0	1	0
63	0	0	1	0	0	0
68	0	0	1	0	0	0
71	0	0	1	0	0	0
82	0	1	2	0	0	0
83	0	0	1	0	0	0
84	0	0	1	0	0	0
96	0	0	1	0	0	0
97	0	1	0	0	0	0
111	0	0	1	0	0	0
117	0	0	1	0	0	0
126	0	0	1	0	0	0
143	0	1	0	0	0	0
148	1	0	0	0	0	0
153	0	1	0	0	0	1
237	1	0	0	0	0	0
360	0	0	0	1	0	0
437	0	1	0	0	0	0
466	0	0	1	0	0	0
486	1	0	0	0	0	0
535	0	0	0	0	1	0
596	0	1	0	0	0	0
647	0	1	0	0	0	0
683	0	1	0	0	0	0
1254	0	0	0	1	0	0
1306	0	1	0	0	0	0
1318	0	0	0	0	1	0
2010	0	0	0	1	0	0

Frequência de ocorrência do num. de falhas por classe de tamanho da falha e mês. A primeira coluna da tabela indica o tamanho da falha.

freq_m <- g %$% 
             table(length, month(sdate)) %>% 
             addmargins(margin = 2)
kable(freq_m, 
      row.names = TRUE, 
      caption = "Tabela 13. Frequência de ocorrência de falhas por classes de tamanho das falhas e mês.  A primeira coluna da tabela indica o tamanho da falha.")

Tabela 13. Frequência de ocorrência de falhas por classes de tamanho das falhas e mês. A primeira coluna da tabela indica o tamanho da falha.
	1	2	3	4	5	6	7	8	9	10	11	12	Sum
1	174	66	92	92	92	74	77	52	57	98	94	140	1108
2	70	44	53	58	38	43	40	27	40	66	47	65	591
3	40	14	33	29	22	26	26	13	17	35	33	27	315
4	18	12	24	14	17	10	10	13	11	25	19	24	197
5	15	10	11	7	10	10	12	8	11	14	17	14	139
6	10	8	7	12	11	10	8	7	10	16	8	16	123
7	12	7	9	12	9	2	7	4	9	5	11	13	100
8	0	7	3	5	4	4	2	2	5	6	5	8	51
9	7	6	7	6	8	4	3	2	3	2	1	12	61
10	2	3	7	3	3	1	1	2	2	8	1	8	41
11	5	3	3	1	3	1	3	1	4	1	3	5	33
12	3	2	0	1	2	1	3	2	0	3	3	3	23
13	3	1	2	4	4	0	1	1	4	3	2	7	32
14	3	1	1	1	3	2	1	1	2	4	3	5	27
15	3	2	2	1	2	1	1	0	5	2	1	1	21
16	1	0	6	0	3	2	1	0	0	1	3	0	17
17	1	1	0	0	2	3	1	0	1	1	0	2	12
18	2	1	0	0	1	0	1	0	1	2	2	4	14
19	1	1	2	0	3	0	0	0	1	3	1	2	14
20	1	2	3	0	0	1	0	2	0	2	0	3	14
21	1	2	1	0	1	2	0	0	0	1	1	3	12
22	0	0	0	1	0	1	1	0	2	1	1	0	7
23	1	0	0	1	2	0	0	1	3	1	2	0	11
24	1	0	0	1	0	0	0	0	0	0	2	1	5
25	1	0	0	0	1	0	0	0	0	2	1	2	7
26	0	1	1	1	0	0	1	0	0	0	2	1	7
27	0	0	0	0	1	0	0	1	1	0	4	2	9
28	0	0	1	0	0	1	0	0	2	0	2	1	7
29	2	0	0	1	0	0	0	0	1	1	1	3	9
30	1	1	0	0	0	0	0	0	0	0	0	2	4
31	0	0	1	1	0	0	0	0	1	0	0	1	4
32	0	1	1	0	1	0	0	1	0	0	1	0	5
33	1	0	2	3	0	0	0	0	1	1	0	0	8
34	1	0	0	0	0	0	0	0	0	0	0	2	3
35	1	0	0	0	0	0	0	0	0	0	1	0	2
36	0	1	0	1	0	0	0	0	0	2	0	0	4
37	0	0	1	0	0	0	0	1	0	1	0	0	3
38	1	0	0	0	1	0	0	0	0	0	0	0	2
39	1	0	0	0	0	0	1	0	0	0	0	0	2
40	0	0	1	0	2	0	0	0	0	2	1	0	6
41	1	1	0	0	0	0	0	0	0	0	1	0	3
42	0	0	1	0	0	0	1	1	1	0	0	0	4
43	0	0	1	0	0	0	0	0	0	0	0	1	2
45	0	0	1	0	0	0	0	0	0	0	1	0	2
46	0	0	0	0	0	0	0	1	0	0	0	0	1
50	0	0	0	0	0	0	0	1	0	0	1	0	2
51	0	0	0	0	0	0	0	0	0	1	0	1	2
52	0	1	0	0	0	0	0	0	1	0	0	0	2
55	0	0	0	0	0	0	0	1	0	0	0	0	1
58	0	0	0	0	0	0	1	0	0	1	0	0	2
61	0	0	0	0	0	0	0	0	0	0	0	1	1
62	0	0	0	0	2	0	0	0	0	0	0	0	2
63	0	0	0	0	0	0	0	1	0	0	0	0	1
68	0	0	0	0	0	0	0	1	0	0	0	0	1
71	0	0	0	0	0	0	0	0	0	0	0	1	1
82	0	0	1	0	0	0	0	0	0	1	1	0	3
83	0	0	1	0	0	0	0	0	0	0	0	0	1
84	1	0	0	0	0	0	0	0	0	0	0	0	1
96	0	0	1	0	0	0	0	0	0	0	0	0	1
97	0	0	1	0	0	0	0	0	0	0	0	0	1
111	0	0	0	0	0	1	0	0	0	0	0	0	1
117	0	0	0	0	0	0	0	0	0	1	0	0	1
126	0	0	0	0	0	0	0	0	0	1	0	0	1
143	0	0	0	1	0	0	0	0	0	0	0	0	1
148	0	0	0	0	0	0	0	0	0	1	0	0	1
153	0	0	1	0	0	0	1	0	0	0	0	0	2
237	0	0	0	0	0	0	0	0	0	0	1	0	1
360	0	1	0	0	0	0	0	0	0	0	0	0	1
437	0	0	0	0	1	0	0	0	0	0	0	0	1
466	0	0	0	0	0	0	0	0	1	0	0	0	1
486	0	0	0	0	0	0	0	0	0	0	0	1	1
535	0	1	0	0	0	0	0	0	0	0	0	0	1
596	0	0	0	0	0	0	0	1	0	0	0	0	1
647	0	0	0	1	0	0	0	0	0	0	0	0	1
683	0	0	0	0	0	0	0	0	1	0	0	0	1
1254	0	0	0	0	1	0	0	0	0	0	0	0	1
1306	1	0	0	0	0	0	0	0	0	0	0	0	1
1318	0	0	0	0	0	0	0	0	0	1	0	0	1
2010	0	0	0	0	0	0	0	0	1	0	0	0	1

Determinação da frequência de ocorrência do total de falhas por mês. A primeira coluna da tabela indica o tamanho da falha.

# frequencia anual
freq_tot_month <- g %$% 
                 table(length, month(sdate)) %>%
                 as.matrix() 
# adicionando margins
freq_tot_month_margins <- freq_tot_month %>% 
  rownames() %>% 
   as.integer() * freq_tot_month %>%
   cbind(., freq_tot_month %>% rowSums()) 
freq_tot_month_margins %<>% colSums() %>% rbind(freq_tot_month_margins, .)
kable(freq_tot_month_margins, 
      row.names = TRUE, 
      caption = "Tabela 14. Total de falhas por classes de tamanho das falhas e mês.  A primeira coluna da tabela indica o tamanho da falha.")

Tabela 14. Total de falhas por classes de tamanho das falhas e mês. A primeira coluna da tabela indica o tamanho da falha.
	1	2	3	4	5	6	7	8	9	10	11	12
1	174	66	92	92	92	74	77	52	57	98	94	140	1108
2	140	88	106	116	76	86	80	54	80	132	94	130	1182
3	120	42	99	87	66	78	78	39	51	105	99	81	945
4	72	48	96	56	68	40	40	52	44	100	76	96	788
5	75	50	55	35	50	50	60	40	55	70	85	70	695
6	60	48	42	72	66	60	48	42	60	96	48	96	738
7	84	49	63	84	63	14	49	28	63	35	77	91	700
8	0	56	24	40	32	32	16	16	40	48	40	64	408
9	63	54	63	54	72	36	27	18	27	18	9	108	549
10	20	30	70	30	30	10	10	20	20	80	10	80	410
11	55	33	33	11	33	11	33	11	44	11	33	55	363
12	36	24	0	12	24	12	36	24	0	36	36	36	276
13	39	13	26	52	52	0	13	13	52	39	26	91	416
14	42	14	14	14	42	28	14	14	28	56	42	70	378
15	45	30	30	15	30	15	15	0	75	30	15	15	315
16	16	0	96	0	48	32	16	0	0	16	48	0	272
17	17	17	0	0	34	51	17	0	17	17	0	34	204
18	36	18	0	0	18	0	18	0	18	36	36	72	252
19	19	19	38	0	57	0	0	0	19	57	19	38	266
20	20	40	60	0	0	20	0	40	0	40	0	60	280
21	21	42	21	0	21	42	0	0	0	21	21	63	252
22	0	0	0	22	0	22	22	0	44	22	22	0	154
23	23	0	0	23	46	0	0	23	69	23	46	0	253
24	24	0	0	24	0	0	0	0	0	0	48	24	120
25	25	0	0	0	25	0	0	0	0	50	25	50	175
26	0	26	26	26	0	0	26	0	0	0	52	26	182
27	0	0	0	0	27	0	0	27	27	0	108	54	243
28	0	0	28	0	0	28	0	0	56	0	56	28	196
29	58	0	0	29	0	0	0	0	29	29	29	87	261
30	30	30	0	0	0	0	0	0	0	0	0	60	120
31	0	0	31	31	0	0	0	0	31	0	0	31	124
32	0	32	32	0	32	0	0	32	0	0	32	0	160
33	33	0	66	99	0	0	0	0	33	33	0	0	264
34	34	0	0	0	0	0	0	0	0	0	0	68	102
35	35	0	0	0	0	0	0	0	0	0	35	0	70
36	0	36	0	36	0	0	0	0	0	72	0	0	144
37	0	0	37	0	0	0	0	37	0	37	0	0	111
38	38	0	0	0	38	0	0	0	0	0	0	0	76
39	39	0	0	0	0	0	39	0	0	0	0	0	78
40	0	0	40	0	80	0	0	0	0	80	40	0	240
41	41	41	0	0	0	0	0	0	0	0	41	0	123
42	0	0	42	0	0	0	42	42	42	0	0	0	168
43	0	0	43	0	0	0	0	0	0	0	0	43	86
45	0	0	45	0	0	0	0	0	0	0	45	0	90
46	0	0	0	0	0	0	0	46	0	0	0	0	46
50	0	0	0	0	0	0	0	50	0	0	50	0	100
51	0	0	0	0	0	0	0	0	0	51	0	51	102
52	0	52	0	0	0	0	0	0	52	0	0	0	104
55	0	0	0	0	0	0	0	55	0	0	0	0	55
58	0	0	0	0	0	0	58	0	0	58	0	0	116
61	0	0	0	0	0	0	0	0	0	0	0	61	61
62	0	0	0	0	124	0	0	0	0	0	0	0	124
63	0	0	0	0	0	0	0	63	0	0	0	0	63
68	0	0	0	0	0	0	0	68	0	0	0	0	68
71	0	0	0	0	0	0	0	0	0	0	0	71	71
82	0	0	82	0	0	0	0	0	0	82	82	0	246
83	0	0	83	0	0	0	0	0	0	0	0	0	83
84	84	0	0	0	0	0	0	0	0	0	0	0	84
96	0	0	96	0	0	0	0	0	0	0	0	0	96
97	0	0	97	0	0	0	0	0	0	0	0	0	97
111	0	0	0	0	0	111	0	0	0	0	0	0	111
117	0	0	0	0	0	0	0	0	0	117	0	0	117
126	0	0	0	0	0	0	0	0	0	126	0	0	126
143	0	0	0	143	0	0	0	0	0	0	0	0	143
148	0	0	0	0	0	0	0	0	0	148	0	0	148
153	0	0	153	0	0	0	153	0	0	0	0	0	306
237	0	0	0	0	0	0	0	0	0	0	237	0	237
360	0	360	0	0	0	0	0	0	0	0	0	0	360
437	0	0	0	0	437	0	0	0	0	0	0	0	437
466	0	0	0	0	0	0	0	0	466	0	0	0	466
486	0	0	0	0	0	0	0	0	0	0	0	486	486
535	0	535	0	0	0	0	0	0	0	0	0	0	535
596	0	0	0	0	0	0	0	596	0	0	0	0	596
647	0	0	0	647	0	0	0	0	0	0	0	0	647
683	0	0	0	0	0	0	0	0	683	0	0	0	683
1254	0	0	0	0	1254	0	0	0	0	0	0	0	1254
1306	1306	0	0	0	0	0	0	0	0	0	0	0	1306
1318	0	0	0	0	0	0	0	0	0	1318	0	0	1318
2010	0	0	0	0	0	0	0	0	2010	0	0	0	2010
.	2924	1893	1929	1850	3037	852	987	1502	4292	3387	1856	2630	27139

8 Escolha das flags de qualidade (QF) dos dados

Variáveis com flag de qualidade rigorosa (fr):

NEE
PAR, Rg

Variáveis com flag de qualidade flexível:

Ustar

Variáveis sem flags de qualidade

lai, swctop, swc1m e swc2m
VPD, Tair

Mantendo somente as flags dos filtros selecionados.

nee_hh_r_qc_final <- nee_hh_r_qc %>% 
  mutate(qc_Ustar_ff = ifelse(Ustar > 2, 1, qc_Ustar_ff)) %>%
  select(date:daytime, qc_NEE_fr, qc_PAR_fr, qc_Rg_fr, qc_Ustar_ff)
glimpse(nee_hh_r_qc_final)

Observations: 87312
Variables:
$ date        (time) 2009-07-08 00:30:00, 2009-07-08 01:00:00, 2009-07...
$ NEE         (dbl) 2.1149520, 4.3855706, 2.1149520, 2.4004219, 1.7394...
$ PAR         (dbl) 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0...
$ Rg          (dbl) 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.000...
$ Tair        (dbl) 17.370, 16.943, 16.870, 16.700, 16.653, 16.537, 16...
$ VPD         (dbl) 4.226812, 3.807304, 3.812228, 3.677811, 3.718069, ...
$ Ustar       (dbl) 0.580, 0.645, 0.618, 0.533, 0.520, 0.504, 0.498, 0...
$ lai         (dbl) 1.5226, 1.5226, 1.5226, 1.5226, 1.5226, 1.5226, 1....
$ swctop      (dbl) 13.6185, 13.6185, 13.6185, 13.6185, 13.6185, 13.61...
$ swc1m       (dbl) 46.76525, 46.76525, 46.76525, 46.76525, 46.76525, ...
$ swc2m       (dbl) 123.2052, 123.2052, 123.2052, 123.2052, 123.2052, ...
$ PotRad      (dbl) 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.000...
$ daytime     (chr) "night", "night", "night", "night", "night", "nigh...
$ qc_NEE_fr   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_PAR_fr   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_Rg_fr    (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ qc_Ustar_ff (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...

Salvando dados com flags de controle de qualidade.

saveRDS(nee_hh_r_qc_final, file = "output/nee_hh_r_qc_flags_pdg.rds")

9 Referências

[1] E. Falge, D. Baldocchi, R. Olson, P. Anthoni, et al. “Gap filling strategies for defensible annual sums of net ecosystem exchange”. In: Agricultural and Forest Meteorology 107.1 (Mar. 2001), pp. 43-69. DOI: 10.1016/s0168-1923(00)00225-2. .

[2] M. Reichstein, E. Falge, D. Baldocchi, D. Papale, et al. “On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm”. In: Global Change Biology 11.9 (2005), pp. 1424-1439. DOI: 10.1111/j.1365-2486.2005.001002.x. .

[3] M. Reichstein and A. M. Moffat. REddyProc: Data processing and plotting utilities of (half-)hourly eddy-covariance measurements. R package version 0.6-0/r9. 2014. .

10 Funções para manipulação dos dados

10.1 `to_oa()`

O pacote openair fornece funções para facilmente gerar gráficos com versatilidade. Para usar essas funções os dados devem ter uma coluna denominada date no formato POSIXt. A função to_oa() faz essa alteração, possibilitando o uso das funções do openair facilmente. Note que a coluna com informação temporal na classe de dados do Pacote REddyProc é também no formato POSIXt porém o nome da variável é DateTime. A função to_oa() faz essa alteração.

10.2 `regularize()`

Os dados podem conter descontinuidades cronológicas (p.ex: uma data faltante nos dados). A função regularize() faz a regularização cronológica (dados com intervalo de tempo regularmente espaçados) da série de dados. A série de dados é expandida de forma que o primeiro dia inicie às 00:30 horas e último dia termine às 00:00 horas. Nesses casos os valores das variáveis sem dados nas datas faltantes ou no período expandido recebem NA.

10.3 `qc_rcr()`

A função qc_rcr() (controle de qualidade baseado na diferença instantânea do NEE com a amplitude do ciclo mediano móvel) foi criada para indicar potenciais dados espúrios das variáveis baseado na comparação entre dois valores:

A amplitude do ciclo diário mediano móvel \(\Delta x_{i}\) com janela de N dias (p.ex.: 60 dias) em torno da meia-hora considerada; e
A diferença \(\delta x_{i}\) do valor da variável de interesse (p.ex.: NEE) entre instante atual \(x_{i}\) e o instante anterior \(x_{i-1}\) .
- Se \(\delta x_{i} > \Delta x_{i}\) o valor daquela meia-hora e marcado como potencial dados espúrio.

A função qc_rcr() foi criada para aplicação do flag de controle de qualidade para dados espúrios. A função possui um parâmetros adicionais:

window.size: tamanho da janela usado para cálculo do ciclo diário médio móvel
mean.function: função para calcular a medida central do valor de cada meia-hora, p.ex.: a média mean; valor default: median (mediana, menos sensível a dados espúrios).
sd.function: função para calcular a medida de dispersão do valor de cada meia-hora, p.ex.: o desvio padrão sd; valor default: mad (desvio absoluto mediano, menos sensível a dados espúrios).
weight.sd.function: fator multiplicativo da amplitude média móvel. A amplitude do ciclo diário médio móvel é determinada pela diferença entre o máximo e o mínimo de \(\bar{x_{i}}\mp w{x_{i}}'\), onde \({x_{i}}'\) é a medida de dispersão (sd.function); \(\bar{x_{i}}\) é a mediana horária da variável e \(w_{i}\) é o peso dado ao resultado da sd.function. Se weight.sd.function = 0 (1) então a amplitude é determinada pela diferença entre os valores máximo e mínimo de \(\bar{x_{i}}\) (\(\bar{x_{i}}\mp w{x_{i}}'\)). Esse parâmetro fornece flexibilidade a filtragem de dados.

A saída dessa função é é um data frame com as colunas: date, var.name (nome da variável, p.ex. NEE, definida pelo parâmetro var.name na chamada da função) e qc_rcr_var.name. Esta última é uma flag que indica se a variável e dado válido (qc_rcr_var.name = 0) ou espúrio (qc_rcr_var.name = 1).

10.4 `qc_hh_quantile()`

O critério do limiar de quantil da meia-hora fornece outra forma para identificar potenciais dados espúrios das variáveis. A função qc_hh_quantile() foi criada para isso e tem a flexibilidade de definir diferentes limiares para os quantis inferior e superior dos dados de meia-hora.

A saída dessa função é é um data frame com as colunas: date, var.name (nome da variável, p.ex. NEE, definida pelo parâmetro var.name na chamada da função) e qc_var.name_qX. Esta última é uma flag que indica se a variável atende (qc_var.name_qX = 0) ou não (qc_var.name_qX = 1) os critérios dos quantis das meia-horas, ou seja o valor da série temporal de 30 min está entre os limiares dos quantis (entre as linhas azul e vermelha nos gráficos abaixo) ou não, respectivamente. Portanto qc_var.name = 1 indica que o valor da variável deve ser filtrado. O sufixo qX no nome da flag (qc_var.name_qX) é resultante do valor do quantil inferior multiplicado por 100. Portanto se o quantil inferior de entrada da função for 0.01, e.g. qc_hh_quantile(X, var.name = "NEE", qtls = c(0.01, 0.98)) a flag de saída será qc_NEE_q1, mesmo que o quantil superior seja um valor que não corresponda ao quantil complementar (0.99).

Testes de filtragem podem ser feitos para fins de comparação usando dois valores de quantil diferentes, por exemplo:

filtragem flexível (ff): com percentis de 0.001 e 0.999 (0.1% dos dados de cada meia-hora, localizados nos extremos da distribuição da variável, serão indicados).
filtragem rigorosa (fr): com percentis de 0.01 e 0.99 (1% dos dados de uma dada meia-hora, localizados nos extremos da distribuição da variável, serão indicados).

Controle de qualidade e levantamento de falhas

Jônatan Tatsch (Departamento de Física - UFSM)

03-06-2015

1 Introdução

2 Pré-configuração: pacotes e scripts necessários

3 Carregando dados

4 Definiçao de períodos diurnos e noturnos

5 Controle de qualidade

5.1 Critérios para indicação de dados espúrios

5.1.1 Aplicação ao `NEE`

5.1.1.1 Filtro baseado na amplitude do ciclo diário

5.1.1.2 Filtro baseado no quantil da meia-hora

5.1.2 Demais variáveis

5.1.3 Série temporal de `NEE` de 30 min com indicação de dados espúrios

5.1.4 Série temporal de `NEE` por hora.

5.1.5 Série temporal de `NEE` médio noturno para os critérios de filtragem

5.1.5.1 Gráfico da Frequência de NEE noturno disponível

5.1.5.2 Gráficos das médias noturnas de NEE por dia.

6 Frequência de ocorrência de dados faltantes e rejeitados

6.1 Dia e Noite

7 Tamanho das falhas

8 Escolha das flags de qualidade (QF) dos dados

9 Referências

10 Funções para manipulação dos dados

10.1 `to_oa()`

10.2 `regularize()`

10.3 `qc_rcr()`

10.4 `qc_hh_quantile()`

Controle de qualidade e levantamento de falhas

Jônatan Tatsch (Departamento de Física - UFSM)

03-06-2015

1 Introdução

2 Pré-configuração: pacotes e scripts necessários

3 Carregando dados

4 Definiçao de períodos diurnos e noturnos

5 Controle de qualidade

5.1 Critérios para indicação de dados espúrios

5.1.1 Aplicação ao NEE

5.1.1.1 Filtro baseado na amplitude do ciclo diário

5.1.1.2 Filtro baseado no quantil da meia-hora

5.1.2 Demais variáveis

5.1.3 Série temporal de NEE de 30 min com indicação de dados espúrios

5.1.4 Série temporal de NEE por hora.

5.1.5 Série temporal de NEE médio noturno para os critérios de filtragem

5.1.5.1 Gráfico da Frequência de NEE noturno disponível

5.1.5.2 Gráficos das médias noturnas de NEE por dia.

6 Frequência de ocorrência de dados faltantes e rejeitados

6.1 Dia e Noite

7 Tamanho das falhas

8 Escolha das flags de qualidade (QF) dos dados

9 Referências

10 Funções para manipulação dos dados

10.1 to_oa()

10.2 regularize()

10.3 qc_rcr()

10.4 qc_hh_quantile()

5.1.1 Aplicação ao `NEE`

5.1.3 Série temporal de `NEE` de 30 min com indicação de dados espúrios

5.1.4 Série temporal de `NEE` por hora.

5.1.5 Série temporal de `NEE` médio noturno para os critérios de filtragem

10.1 `to_oa()`

10.2 `regularize()`

10.3 `qc_rcr()`

10.4 `qc_hh_quantile()`