Preparação de dados da Torre-Micro LAI/FPAR
Roilan Hernández Valdés
Dados disponiveis são, LAI, FPAR e Qualidade dos mesmos. Essas informações foram obtidas para o diminio da bacia Pauliceia (reprojetadas e recortadas para o domínio). Em formato raster:
# Lendo data com
lai.fpar.qc.rasters <- readRDS("./rastersLAIFPARQC.rds")
print(lai.fpar.qc.rasters)## [[1]]
## class : RasterStack
## dimensions : 14, 12, 168, 1141 (nrow, ncol, ncell, nlayers)
## resolution : 500, 500 (x, y)
## extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
## names : LAI_2002.08.05, LAI_2002.08.09, LAI_2002.08.13, LAI_2002.08.17, LAI_2002.08.21, LAI_2002.08.25, LAI_2002.08.29, LAI_2002.09.02, LAI_2002.09.06, LAI_2002.09.10, LAI_2002.09.14, LAI_2002.09.18, LAI_2002.09.22, LAI_2002.09.26, LAI_2002.09.30, ...
## min values : 0.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ...
## max values : 255.000000, 67.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, 255.000000, ...
##
##
## [[2]]
## class : RasterStack
## dimensions : 14, 12, 168, 1141 (nrow, ncol, ncell, nlayers)
## resolution : 500, 500 (x, y)
## extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
## names : LAI_qc_2002.08.05, LAI_qc_2002.08.09, LAI_qc_2002.08.13, LAI_qc_2002.08.17, LAI_qc_2002.08.21, LAI_qc_2002.08.25, LAI_qc_2002.08.29, LAI_qc_2002.09.02, LAI_qc_2002.09.06, LAI_qc_2002.09.10, LAI_qc_2002.09.14, LAI_qc_2002.09.18, LAI_qc_2002.09.22, LAI_qc_2002.09.26, LAI_qc_2002.09.30, ...
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## max values : 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, ...
##
##
## [[3]]
## class : RasterStack
## dimensions : 14, 12, 168, 1141 (nrow, ncol, ncell, nlayers)
## resolution : 500, 500 (x, y)
## extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
## names : FPAR_2002.08.05, FPAR_2002.08.09, FPAR_2002.08.13, FPAR_2002.08.17, FPAR_2002.08.21, FPAR_2002.08.25, FPAR_2002.08.29, FPAR_2002.09.02, FPAR_2002.09.06, FPAR_2002.09.10, FPAR_2002.09.14, FPAR_2002.09.18, FPAR_2002.09.22, FPAR_2002.09.26, FPAR_2002.09.30, ...
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## max values : 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, ...Informação de qualidade dos dados de FPAR e LAI, tem a seguente estrutura:
#
# FparLai_QC 5 BITFIELDS IN 8 BITWORD
# MODLAND_QC START 0 END 0 VALIDS 2
# MODLAND_QC 0 = Good Quality (main algorithm with or without saturation)
# MODLAND_QC 1 = Other Quality (back-up algorithm or fill value)
# SENSOR START 1 END 1 VALIDS 2
# SENSOR 0 = Terra
# SENSOR 1 = Aqua
# DEADDETECTOR START 2 END 2 VALIDS 2
# DEADDETECTOR 0 = Detectors apparently fine for up to 50% of channels 1,2
# DEADDETECTOR 1 = Dead detectors caused >50% adjacent detector retrieval
# CLOUDSTATE START 3 END 4 VALIDS 4 (this inherited from Aggregate_QC bits {0,1} cloud state)
# CLOUDSTATE 00 = 0 Significant clouds NOT present (clear)
# CLOUDSTATE 01 = 1 Significant clouds WERE present
# CLOUDSTATE 10 = 2 Mixed cloud present on pixel
# CLOUDSTATE 11 = 3 Cloud state not defined,assumed clear
# SCF_QC START 5 END 7 VALIDS 5
# SCF_QC 000=0 Main (RT) algorithm used, best result possible (no saturation)
# SCF_QC 001=1 Main (RT) algorithm used, saturation occured. Good, very usable.
# SCF_QC 010=2 Main algorithm failed due to bad geometry, empirical algorithm used
# SCF_QC 011=3 Main algorithm faild due to problems other than geometry, empirical algorithm used
# SCF_QC 100=4 Pixel not produced at all, value coudn't be retrieved (possible reasons: bad L1B data,
# unusable MODAGAGG data) # Vetor de datas extraído do nomes dos layers do arquivos
# LAI
dates <- as.Date(substr(names(lai.fpar.qc.rasters[[1]]),5,14), format="%Y.%m.%d")
# Correção aos dados de LAI brutos
lai.list <- llply(1:nlayers(lai.fpar.qc.rasters),.progress = "text",
function(x){ # x <- 43
# Eliminando os dados fora do intervalo válido de LAI
lai.rasters[[x]][lai.rasters[[x]] > 100] <- NA
# Preenchendo as células com NA, a matriz de convolução ..
# .. tem que ser de 5x5 para conseguir preencher todas ..
# .. as células com NA, é caso particular.
lai.rasters[[x]] <- focal(lai.rasters[[x]],
pad=TRUE,
w = matrix(1/25,5,5),
NAonly=TRUE,
na.rm=TRUE)
# Aplicando correção dos valores de LAI, para..
# .. levar ao intervalo de 0,10. (o real é de (0,7))
gain(lai.rasters[[x]]) <- 0.1
return(lai.rasters[[x]])
}
)
# Criando data.frame com série temporal ao longo das colunas..
# .. e em cada linha um ponto da grade
lai.val <- llply(lai.list,.progress = "text",
function(i){as.numeric(getValues(i))}
) %>% as.data.frame %>% setNames(dates)
# Criando data.frame com série temporal ao longo das colunas..
# .. e em cada linha um ponto da grade
qc.list <- unstack(lai.qc.rasters)
qc.val <- llply(qc.list,.progress = "text",
function(i){getValues(i)}
) %>% as.data.frame %>% setNames(dates) # Determinando continuidade da série
dates_seq <- seq.Date(from = dates[1],to = dates[length(dates)], by = "4 day")
splitdates <- split(dates, years(dates))
# dates list sem datas faltantes com intervalo
# dos arquivos de LAI
dates_l <- lapply(splitdates, function(x) {
merge(data.frame(date = seq(x[1], x[length(x)], by="4 day")),
data.frame(date=x, difdate=c(4,diff(x))),
all = T)
} ) %>% bind_rows
# Datas faltantes??
date_falt <- dates_l$date[!(dates_l$date %in% dates)]
date_falt## [1] "2003-12-19" # Interpolação linear para datas faltantes.
# Preenchimento do
# >> lai.val nas datas faltantes com interpolação linear
# >> qc.val com as flags do dia anterior (não é possivel interpolar)
for(i in 1:length(date_falt)){
a <- which(dates_l$date == date_falt[i])
b <- data.frame((lai.val[,a-1]+lai.val[,a])/2) %>% setNames(date_falt[i])
lai.val <- cbind(lai.val[,1:(a-1)],b,lai.val[,a:ncol(lai.val)])
qc.val <- cbind(qc.val[,1:(a-1)],qc.val[,a],qc.val[,a:ncol(qc.val)])
}# Por requerimentos do TIMESAT, os dados devem fechar anos completos
# Faremos então o completamento das séries de LAI e FPAR com o ciclo médio do LAI
# Definindo sequência de datas a ser preenchidas
comp.dates <- seq.Date(from = as.Date("2002-01-01"),to = dates[1], by = "4 day")
# Criando data frame com completamento
complemento <- data.frame(matrix(NA,ncol = length(comp.dates)-1,nrow = length(lai.val[,1])))
# Nomeando novas colunas com a data
names(complemento) <- comp.dates[-length(comp.dates)]
# Juntando data frames de completamento e dados originais
# para LAI...
lai.val.comp <- cbind(complemento,lai.val)
# para QC.
qc.val.comp <- cbind(complemento,qc.val) # dados uma serie temporal (padronizada, i.e. com 1196 ptos no caso pauliceia,
# determina o ciclo anual medio e preenche ptos faltantes da serie de entrada
# com o valor medio do ciclo anual)
# (FUNÇÃO criada por TATSCH, 20xx)
getcycleAvgTS <- function(x, each8day) {
viclima <- round(rep(c(t(unlist(lapply(split(t(x),each8day), mean, na.rm=T)))),13), 2)
data2fill <- viclima[is.na(x)]
return(data2fill)
} # vetor indicando ciclo das series
each4day <- rep(1:92, 13)
# preenchendo pontos das series correspondentes as novas datas
# inseridas para padronizacao da serie
for(i in 1:nrow(lai.val.comp)) {
# cat(i, " de "," ",nrow(lai.val.comp), "\n")
lai.val.comp[i, is.na(lai.val.comp[i,])] <- getcycleAvgTS(x = lai.val.comp[i, ], each8day = each4day)
} # Função para converão de úmero decimal a binario de 8 bits
to8bitword <- function(x=157,nbits = 8){
sapply(x,function(i){
paste(rev(as.integer(intToBits(i))), collapse="") %>% substr((nbits-1),32)
}) %>% as.vector
# return()
}
# Conversão dos dados de qualidade a peso para utilização do TIMESAT.
for(i in 1:nrow(lai.val.comp)) { # i=1
##cat(i, " de "," ",nrow(lai.val.comp), "\n")
# Os valores de NA no qc.val.comp correspondem ao ciclo médio anual,
# Podemos preencher com 0, a 8bitsword == "00000000"
qc.val.comp[i, is.na(qc.val.comp[i,])] <- 0
# Criando vetor de peso 0 para todos os pontos inicialmente
qc.peso <- rep(0,times=length(qc.val.comp[i,]))
# Peso =1, onde o 8bitword do valor de QC nas posições 4 e 5 é: "00" ou "11"
# CLOUDSTATE 00 = 0 Significant clouds NOT present (clear)
# CLOUDSTATE 01 = 1 Significant clouds WERE present
# CLOUDSTATE 10 = 2 Mixed cloud present on pixel
# CLOUDSTATE 11 = 3 Cloud state not defined,assumed clear
qc.peso[which(substr(to8bitword(qc.val.comp[i,]),4,5) %in% c("00","11"))] <- 1
qc.peso[which(substr(to8bitword(qc.val.comp[i,]),6,8) %in% c("010","011","100") & qc.peso == 1)] <- 0.5
# qc.peso[which(substr(to8bitword(qc.val.comp[i,]),3,3) == "1")] <- 0
qc.val.comp[i,] <- qc.peso
} # Escrevendo os arquivos no formato TIMESAT
output_dir <- "~/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran/"
if(!dir.exists(output_dir))dir.create(output_dir)
# Cabeçalhos >>
write.fwf(file=paste0(output_dir,"LAISat.txt"),x = data.frame(NA,13, 92, 168),
rownames = F,colnames = F,na = "")
write.fwf(file=paste0(output_dir,"QCSat.txt"),x = data.frame(NA,13, 92, 168),
rownames = F,colnames = F,na = "")
# Dados >>
write.fwf(file=paste0(output_dir,"LAISat.txt"),
x= data.frame(dummy=NA,round(lai.val.comp,2)),
rownames = F,colnames = F,na = "",append = TRUE)
write.fwf(file=paste0(output_dir,"QCSat.txt"),
x= data.frame(dummy=NA,qc.val.comp),
rownames = F,colnames = F,na = "",append = TRUE)# Função para escrever configuração de parámetros do TIMESAT
writeTSS <- function(dir = "~/Dissertação/MODIS/TIMESAT_in/",
set = "input_lai.txt",
file = "LAISat.txt",
range = c(0.0,7.0),
mask = "QCSat.txt",
qcvetor = c(0,0,0,0.5,0.5,0.5,1,1,1),
A = 0,
S = 2.0,
season = 0.5,
fit.steps = 3,
adapt = 2,
methods = c(1,1,1),
windw = c(4,5,6),
fseasampl = 30,
jobname="pauliceia"
){
sink(paste0(dir,set))
cat(file,"\n") # LAISat.txt
cat(range,"\n") # 0.0000 7.0000
cat(mask,"\n") # NONE
cat(qcvetor,"\n") # 0,0,0,1,1,1,2,100,0,101,256,0
cat(A,"\n") # 0.0000
cat(S,"\n") # 2.0000
cat(season,"\n") # 0.5000
cat(fit.steps,"\n") # 3
cat(adapt,"\n") # 2.0000
cat(methods,"\n") # 1 1 1
cat(windw,"\n") # 4 5 6
cat(fseasampl,"\n") # 20.0000
cat(jobname,"\n") # pauliceia
sink()
} # Teste para os parâmetros de entrada ao TIMESAT
list(S = seq(1,5,0.5),
# season = seq(0,1,0.1),
# fit.steps = c(1:3),
adapt = c(1:10),
# methods = c(1,1,1),
windw = c(4,8,12), #c(1,1.5,2.5)
fseasampl = seq(10,60,10)
) -> param.setllply(1:length(param.set), function(i){
# i = 4
sapply(1:length(param.set[[i]]), function(j){
# j = 4
if(names(param.set)[i] == "S") writeTSS(dir = "~/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran/",
S = param.set[["S"]][j])
if(names(param.set)[i] == "adapt")writeTSS(dir = "~/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran/",
adapt = param.set[["adapt"]][j])
if(names(param.set)[i] == "windw")writeTSS(dir = "~/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran/",
windw = c(1,2,3)*param.set[["windw"]][j])
if(names(param.set)[i] == "fseasampl")writeTSS(dir = "~/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran/",
fseasampl = param.set[["fseasampl"]][j])
cat("Running with ",names(param.set)[i], "=", param.set[[i]][j], "\n" )
system(command = paste0("cd /home/roilan/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran",
";./a.out < input_lai.txt"),intern =FALSE,ignore.stdout = TRUE)
output_dir_out <- "~/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran/"
lai.fit.SG <- read.table(file = paste0(output_dir_out,"fort.917"),header = F)
as.numeric(lai.fit.SG[30,])
}) %>% as.data.frame()
}) -> a
names(a) <- names(param.set)’ Single spikes are detected by a comparison with median filtered values and with closest neighbors. If the distance is greater than S*ystd, where ystd is the standard deviation for data values, we have a spike. S = 2 is the normal value.’
Menor valor de S reduz os picos nos dados.
## Warning in RColorBrewer::brewer.pal(length(param.set[[var]]), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colorsGive the strength of the adaptation. Value should be in the range [1,10] where 2 is the normal value.
Maior ADAPT produz um maior ajuste aos máximos e afeta os mínimos que acontecem no inverno.
## Warning in RColorBrewer::brewer.pal(length(param.set[[var]]), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colorsSavitzky-Golay window size for each of the fitting steps NOTA: no caso é definido a janela para o primeiro passo, as janelas para os seguentes casos são definidos pela multiplicação com o vetor
1, 1.5, 2(ex., dada janela inicial 4 a sequência de janelas para os três passos seria: 4, 6, 8)
Janela de menor dimensão não elimina o ruido dos dados, a de maior dimensão não representa os mínimos do inverno. A proposta de [8,12,16] pode ser um boa escolha.
The time for the season start (end) is defined as the time for which the sensor data value ,measured from the base level, has increased (decreased) to X % of the seasonal amplitude. X = 20 is the normal value.
A amplitude sazonal, não mostrou sensibilidade.
Filtragem proposto para a célula da Torre Micrometeorologica
# Criando o arquivo de configurações para o TIMESAT
writeTSS(dir = "~/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran/",
file = "LAISat.txt", # Arquivo com o valor de LAI
range = c(0.0,7.0), # Intervalo válido do LAI
mask = "QCSat.txt", # Arquivo com flags de qualidade
A = 0, # Garanta o processamento de toda a série
methods = c(1,1,1), # Switch dos métodos de suavização
fit.steps = 3, # Quantidade de passos no método Savitzky-Golay
qcvetor = c(0,0,0,0.5,0.5,0.5,1,1,1), # intervalos e peso para cada um...
# [a1,b1,w1....a3,b3,w3] se a1 < qc < b1 então peso = w1 ....
S =1, # Critério para determinação dos picos
season = 0.5, # Parâmetro para determinar número de estaçoes no ano
adapt = 5, # Nível de forçamento ao envelope máximo da série
windw = c(1,1.5,2)*8, # Janela temporal para cada passo no método Savitzky-Golay
fseasampl = 20 # Porcentagem da amplitude sazonal de LAI para determinar o ...
# ... o momento de transição de estação.
)
system(command = paste0("cd /home/roilan/Dropbox/Dissertacao/MODIS/timesat2_3/timesat_fortran/",
";./a.out < input_lai.txt"),intern = TRUE)