Nothing
#' curvefits by local model functions of TIMESAT
#'
#' @description Local model functions `f_L(t)`, `f_C(t)` and `f_R(t)`
#' describe the VI variation in intervals around the left minima, the central
#' maxima and the right minima.
#'
#' Local model function are merged into global model function via [merge_LocalModels()]
#' and Per J\"onsson et al. (2004; their Eq. 12),
#' where cut-off function sharply drop from 1 to 0 in small intervals around
#' `(t_L + t_C)/2` and `(t_C + t_R)/2`.
#'
#' \deqn{
#' F(t)= \begin{cases}
#' \alpha(t) f_{L}(t)+[1-\alpha(t)] f_{C}(t), t_{L}<t<t_{C} \\
#' \beta(t) f_{C}(t)+[1-\beta(t)] f_{R}(t), t_{C}<t<t_{R}\end{cases}
#' }
#'
#' @inheritParams curvefits
#' @param fits List objects returned by [curvefits_LocalModel()] (not [curvefits()]).
#'
#' @inheritSection curvefits options for fitting
#'
#' @references
#' 1. Per J\"onsson, P., Eklundh, L., 2004. TIMESAT - A program for analyzing
#' time-series of satellite sensor data. Comput. Geosci. 30, 833-845.
#' \doi{10.1016/j.cageo.2004.05.006}.
#'
#' @seealso [curvefits()]
#' @example R/examples/ex-curvefits_LocalModel.R
#' @export
curvefits_LocalModel <- function(
INPUT, brks,
# methods, wFUN,
options = list(),
# iters = 2, wmin = 0.1,
# nextend = 2, maxExtendMonth = 2, minExtendMonth = 1,
# minT = 0,
# minPercValid = 0,
# use.rough = FALSE,
# use.y0 = TRUE,
...)
{
if (all(is.na(INPUT$y))) return(NULL)
set_options(fitting = options)
opt = .options$fitting
QC_flag <- INPUT$QC_flag
nptperyear <- INPUT$nptperyear
t <- INPUT$t
years <- year(t)
n <- length(t)
# doys <- as.numeric(t) # days from origin
doys <- as.numeric(difftime(t, date.origin, units = "day")) # + 1
# Tn for background module
w <- w0 <- INPUT$w
y0 <- if (opt$use.y0) INPUT$y0 else INPUT$y
# Tn <- INPUT$Tn # if has no Tn, NULL will be return
# has_Tn <- ifelse(is_empty(Tn), FALSE, TRUE)
# possible snow or cloud, replaced with Whittaker smoothing.
I_all <- match(brks$fit$t, t) %>% rm_empty()
if (opt$use.rough) {
# if the range of t is smaller than `fit$t`
INPUT$y[I_all] <- last(brks$fit)
} else {
I_fix <- which(w[I_all] == opt$wmin)
I_y <- I_all[I_fix]
INPUT$y[I_y] <- last(brks$fit)[I_fix]
}
# use the weights of last iteration of rough fitting
# w[I_all] <- brks$fit %>% {.[, contain(., "witer"), with = F]} %>% last()
# w[I_fix] <- wmin + 0.1 # exert the function of whitaker smoother
# growing season division
di <- data.table( beg = getDateId_before(brks$dt$beg, t),
peak = getDateId_before(brks$dt$peak, t),
end = getDateId_after(brks$dt$end, t)) #%>% na.omit()
# TODO: fix the bug when seasons are not continuous
s = c(di$beg, di$peak, di$end %>% last()) %>% sort()
nlocal = length(s) - 2
width_ylu = nptperyear*2
y <- INPUT$y
fits <- vector(nrow(di), mode = "list")
curvefit_gs <- function(I_beg, I_end, type = 1L) {
I <- I_beg:I_end
# I_extend <- get_extentI(w0, I_beg, I_end, nptperyear)
# In TIMESAT global model function, `nextent` is not used.
I_extend <- I
## 2. input data
ti <- doys[I_extend]
yi <- INPUT$y[I_extend]
wi <- w[I_extend]
# yi_good <- yi[w0[I_extend] > wmin]
## update ylu in a three year moving window
ylu <- get_ylu(y, years, w0, width_ylu, I_beg:I_end, Imedian = TRUE, opt$wmin)
ylu <- merge_ylu(INPUT$ylu, ylu)
# yi[yi < ylu[1]] <- ylu[1] # update y value
# beginI <- ifelse(i == 1, 1, 2) # make sure no overlap
beginI = 1
tout <- doys[I] %>% {.[beginI]:last(.)} # make sure return the same length result.
fit <- curvefit(yi, ti, tout, nptperyear = nptperyear,
w = wi, ylu = ylu,
iters = opt$iters, methods = opt$methods, wFUN = opt$wFUN,
type = type,
...)
# if too much missing values
if (sum(wi > pmax(opt$wmin+0.1, 0.2))/length(wi) < opt$minPercValid){
fit$model %<>% map(function(x){
x$zs %<>% map(~.x*NA) # list()
return(x)
})
}
return(fit)
}
fits = list()
for (i in 1:nlocal) {
# odd: t2t (1, 3, 5, ...)
# even: p2k (2, 4, 6)
# typeI = ifelse(i %% 2 == 1, "t2t", "p2p")
typeI = ifelse(i %% 2 == 1, 1L, -1L)
if (opt$verbose) fprintf(" [curvefits] running %d ... \n", i)
I = s[i]:s[i+2]
data = list(t = doys[I], y = y0[I], QC_flag = QC_flag[I]) %>% as.data.table()
fit = curvefit_gs(s[i], s[i+2], type = typeI)
fit$data = data
fit$type = typeI
fits[[i]] = fit
}
names(fits)[seq(1, nlocal, 2)] = brks$dt$flag
return(fits)
# L1:curve fitting method, L2:yearly flag
# return(list(tout = t[first(di$beg):last(di$end)], # dates for OUTPUT curve fitting VI
# fits = fits))
}
#' cutoff
#'
#' @param n1,n2 peak and trough (or trough and peak)
#'
#' @return A numeric vector, in small intervals around `(tL + tC)/2` and `(tC+tR)/2`,
#' respectively, smoothly drop from 1 to 0.
#' @references
#' 1. Per J\"onsson, P., Eklundh, L., 2004. TIMESAT - A program for analyzing
#' time-series of satellite sensor data. Comput. Geosci. 30, 833-845.
#' https://doi.org/10.1016/j.cageo.2004.05.006.
#' @keywords internal
cutoff <- function(n1, n2, k = n1:n2) {
ndiff = n2 - n1;
nmid = (n1 + n2) / 2;
# k = n1:n2;
1 - (atan(10 * (k - nmid) / ndiff) - atan(10 * (n1 - nmid) / ndiff)) /
(atan(10 * (n2 - nmid) / ndiff) - atan(10 * (n1 - nmid) / ndiff))
}
#' @rdname curvefits_LocalModel
#' @export
merge_LocalModels <- function(fits) {
# begin from `t2t`; end at `t2t`
nlocal = length(fits)
N_gs = (nlocal + 1) / 2
fits2 = list()
for (i in 1:N_gs) {
# i_c = c(di$beg[i], di$end[i])
# i_l = c(di$peak[i-1], di$peak[i])
# i_r = c(di$peak[i], di$peak[i+1])
k = (i - 1) * 2 + 1
fit_c = fits[[k]]
fit_l = if(k <= 1) NULL else fits[[k-1]]
fit_r = if(k >= N_gs) NULL else fits[[k+1]]
fit = .merge_LocalModel(fit_c, fit_l, fit_r)
fits2[[i]] <- fit
}
names(fits2) = names(fits)[seq(1, nlocal, 2)]
fits2
}
# ' Merge Local Model Functions
# ' @references
# ' TIMESAT Mannual v3.3 (2018), Eq. 10
# ' @rdname curvefits_LocalModel
# ' @export
.merge_LocalModel <- function(fit_c, fit_l = NULL, fit_r = NULL) {
if (is.null(fit_l) && is.null(fit_r)) return(fit_c)
t_c = fit_c$tout
t_l = fit_l$tout
t_r = fit_r$tout
info_l = match2(t_l, t_c)
info_r = match2(t_r, t_c)
a = if (!is.null(fit_l)) cutoff(first(info_l$I_x), last(info_l$I_x)) else 0
b = if (!is.null(fit_r)) cutoff(first(info_r$I_x), last(info_r$I_x)) else 1
# for z_c
ind_l = info_l$I_y
ind_r = info_r$I_y
if (is.null(fit_l) && !is.null(fit_r)) ind_l = -(ind_r)
if (is.null(fit_r) && !is.null(fit_l)) ind_r = -(ind_l)
nmeth = fit_c$model %>% length()
niter = fit_c$model[[1]]$zs %>% length()
iterator = set_names(1:niter, paste0("iter", 1:niter))
for(i in 1:nmeth) {
zs = map(iterator, function(j){
# for (j in 1:niter) {
z_c = fit_c$model[[i]]$zs[[j]]
z_l = if (is.null(fit_l)) 0 else fit_l$model[[i]]$zs[[j]][info_l$I_x]
z_r = if (is.null(fit_r)) 0 else fit_r$model[[i]]$zs[[j]][info_r$I_x]
z_cl = z_c[ind_l]
z_cr = z_c[ind_r]
# TODO: 目前算法可能会导致zc变短,L84 `beginI = 1`可消除此错误
ans = c(a * z_l + (1 - a) * z_cl,
b * z_cr + (1 - b) * z_r)
# plot(z_c, col = "grey", type = "l", lwd = 1, ylim = c(0, 0.6))
# lines(ind_l, z_l)
# lines(ind_r, z_r)
# lines(ans, col = "red")
ans
})
fit_c$model[[i]]$zs = zs
}
fit_c
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.