Skip to content

Commit 101b557

Browse files
MrJGaoMrJGao
committed
Add greendown_aware for autumn phenometrics retrieval
- Remove unused package dependencies - Account for the summer greendown issue when using thresholds to retreive phenometrics
1 parent bd76800 commit 101b557

4 files changed

Lines changed: 105 additions & 47 deletions

File tree

‎DESCRIPTION‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,6 @@ RoxygenNote: 7.3.1
1717
Imports:
1818
coda,
1919
data.table,
20-
geojsonio,
21-
geojsonR,
22-
getPass,
2320
httr,
2421
jsonlite,
2522
lubridate,

‎R/hlp_funs.R‎

Lines changed: 80 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -279,6 +279,25 @@ CalPhenoParam <- function(
279279
}
280280

281281

282+
#' Find point to line distance
283+
#'
284+
#' @param x1,y1,x2,y2 The coordinates of the two points that define the line.
285+
#' @param x0,y0 The point coordinates to calculate the distance.
286+
#'
287+
#' @return Distance from the point to the line.
288+
#' @noRd
289+
CalPointToLineDistance <- function(x1, y1, x2, y2, x0, y0) {
290+
m <- (y2 - y1) / (x2 - x1)
291+
A <- - m
292+
B <- 1
293+
C <- m * x1 - y1
294+
295+
distance <- abs(A * x0 + B * y0 + C) / sqrt(A^2 + B^2)
296+
297+
return(distance)
298+
}
299+
300+
282301
#' Calculate phenometrics using the threshold-based method
283302
#'
284303
#' @param p_li A list containing the model parameters.
@@ -295,16 +314,20 @@ CalPhenoParam <- function(
295314
#' The default level is 0.9, generating `90%` credible intervals. The end
296315
#' points of these intervals define the upper and lower bounds for the estimated
297316
#' phenometrics.
317+
#' @param greendown_aware Default is `FALSE`. If `TRUE`, Senescence will be
318+
#' retrieved as the end of summer greendown date, and MidGreendown as the
319+
#' transition point of the first derivative of the autumn EVI2 curve.
298320
#'
299-
#' #' @return An "BlspFit" object filled with retrieved phenology and parameters.
321+
#' @return An "BlspFit" object filled with retrieved phenology and parameters.
300322
#'
301323
#' @import data.table
302324
#' @noRd
303325
CalPhenoThresh <- function(
304326
p_li, mod,
305327
years, numYears,
306328
date_vec, vi_vec, weights_vec,
307-
cred_int_level
329+
cred_int_level,
330+
greendown_aware = FALSE
308331
) {
309332
# Format MCD12Q2-like phenometrics
310333
blsp_fit <- EmptyBlspOutput(
@@ -332,14 +355,7 @@ CalPhenoThresh <- function(
332355
# Spring amp
333356
spring_min <- min(p[1:peakdate])
334357
spring_amp <- peak - spring_min
335-
# Autumn amp
336-
autumn_min <- min(p[peakdate:length(p)])
337-
autumn_amp <- peak - autumn_min
338-
339-
if (spring_amp < 0.2 | autumn_amp < 0.2) {
340-
return(rep(NA, 7))
341-
}
342-
358+
343359
gup <- which(
344360
p[1:peakdate] > (spring_amp * 0.15 + min(p[1:peakdate]))
345361
)[1]
@@ -349,16 +365,60 @@ CalPhenoThresh <- function(
349365
mat <- which(
350366
p[1:peakdate] > (spring_amp * 0.90 + min(p[1:peakdate]))
351367
)[1]
352-
353-
sens <- which(
354-
p[peakdate:length(p)] < autumn_amp * 0.90 + autumn_min
355-
)[1] + peakdate
356-
midgdown <- which(
357-
p[peakdate:length(p)] < autumn_amp * 0.5 + autumn_min
358-
)[1] + peakdate
359-
dorm <- which(
360-
p[peakdate:length(p)] < autumn_amp * 0.15 + autumn_min
361-
)[1] + peakdate
368+
369+
if (isTRUE(greendown_aware)) {
370+
# If greendown aware, autumn amp is calculated as the EVI2
371+
# between the EVI2 of the end of greendown to the EVI2 of
372+
# dormancy
373+
374+
# Draw a line between peak and min, find the point w/ the max
375+
# distance to the line, which is senescence
376+
autumn_p <- p[peakdate:length(p)]
377+
autumn_min <- min(autumn_p)
378+
autumn_amp <- peak - autumn_min
379+
autumn_min_doy <- which(
380+
autumn_p < autumn_amp * 0.1 + autumn_min
381+
)[1] + peakdate
382+
383+
dis <- sapply(peakdate + seq_along(autumn_p) - 1, function(x0) {
384+
CalPointToLineDistance(
385+
peakdate, peak,
386+
autumn_min_doy, autumn_min,
387+
x0, autumn_p[x0 - peakdate + 1]
388+
)
389+
})
390+
sens <- peakdate + which(diff(sign(diff(dis))) == -2) + 2
391+
392+
# Now, the autumn amplitude is from sens to min
393+
autumn_amp <- p[sens] - autumn_min
394+
395+
# MidGreendown is from the first derivative
396+
midgdown <- which.min(
397+
diff(autumn_p[(sens - peakdate + 1):length(autumn_p)])
398+
) + sens
399+
dorm <- which(
400+
autumn_p[(sens - peakdate + 1):length(autumn_p)] <
401+
autumn_amp * 0.15 + autumn_min
402+
)[1] + sens
403+
} else {
404+
# Autumn amp
405+
autumn_min <- min(p[peakdate:length(p)])
406+
autumn_amp <- peak - autumn_min
407+
408+
if (spring_amp < 0.2 | autumn_amp < 0.2) {
409+
return(rep(NA, 7))
410+
}
411+
412+
sens <- which(
413+
p[peakdate:length(p)] < autumn_amp * 0.90 + autumn_min
414+
)[1] + peakdate
415+
midgdown <- which(
416+
p[peakdate:length(p)] < autumn_amp * 0.5 + autumn_min
417+
)[1] + peakdate
418+
dorm <- which(
419+
p[peakdate:length(p)] < autumn_amp * 0.15 + autumn_min
420+
)[1] + peakdate
421+
}
362422

363423
return(c(gup, midgup, mat, peakdate, sens, midgdown, dorm))
364424
})

‎R/mod_fit.R‎

Lines changed: 25 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -101,31 +101,34 @@ BLSPFitted <- function(blsp_fit, asCI = FALSE) {
101101
#' to the supplied data (can be sparse), and returns phenometrics for the
102102
#' entire time frame. For further explanation, please see the vignette.
103103
#'
104-
#' @param date_vec The date vector, be sure to convert the vector to "Date"
105-
#' format or use "yyyy-mm-dd" format string.
106-
#' @param vi_vec The vegetation index vector.
104+
#' @param date_vec The date vector, be sure to convert the vector to "Date"
105+
#' format or use "yyyy-mm-dd" format string.
106+
#' @param vi_vec The vegetation index vector.
107107
#' @param weights_vec A numeric vector of same length as vi_vec specifying the
108-
#' weights for the supplied observations. Must be between 0 and 1, inclusive.
109-
#' @param model A string indicating the model name. For now, only support
110-
#' "dblog7" and "dblog6" for the 7- and 6-parameter double-logistic functions.
108+
#' weights for the supplied observations. Must be between 0 and 1, inclusive.
109+
#' @param model A string indicating the model name. For now, only support
110+
#' "dblog7" and "dblog6" for the 7- and 6-parameter double-logistic functions.
111111
#' @param init_values Initial values for MCMC sampling. By default, it is
112112
#' assgined `NULL`. It could also be an object returned from the `FitAvgModel()`
113113
#' function that fits an averaged model or a numeric vector provided by the user.
114-
#' @param start_yr The start year of the result. Default is NULL, which means
115-
#' determined by data.
116-
#' @param end_yr The end year of the result. Default is NULL, which means
117-
#' determined by data.
118-
#' @param cred_int_level A scalar value from 0 to 1 (exclusive) that specifies
119-
#' the level for equal-tailed credible intervals of the estimated phenometrics.
120-
#' The default level is 0.9, generating `90%` credible intervals. The end
121-
#' points of these intervals define the upper and lower bounds for the estimated
122-
#' phenometrics.
123-
#' @param opt An option list that contains additional configurations. For now,
124-
#' only support `list(method = "threshold")` to indicate that use the
125-
#' threshold-based method to retrive phenometrics instead of the default SOS and
126-
#' EOS. The threshold-based method will produce 7 phenometrics including
114+
#' @param start_yr The start year of the result. Default is NULL, which means
115+
#' determined by data.
116+
#' @param end_yr The end year of the result. Default is NULL, which means
117+
#' determined by data.
118+
#' @param cred_int_level A scalar value from 0 to 1 (exclusive) that specifies
119+
#' the level for equal-tailed credible intervals of the estimated phenometrics.
120+
#' The default level is 0.9, generating `90%` credible intervals. The end points
121+
#' of these intervals define the upper and lower bounds for the estimated
122+
#' phenometrics.
123+
#' @param opt An option list that contains additional configurations. It
124+
#' supports `list(method = "threshold")` to indicate that use the
125+
#' threshold-based method to retrive phenometrics instead of the default SOS
126+
#' and EOS. The threshold-based method will produce 7 phenometrics including
127127
#' Greenup, MidGreenup, Maturity, Peak, Senescence, MidGreendown, and Dormancy
128-
#' using VI amplitude thresholds of 15%, 50%, 90%, and 100%, respectively.
128+
#' using VI amplitude thresholds of 15%, 50%, 90%, and 100%, respectively. If
129+
#' including `greendown_aware = TRUE`, Senescence will be retrieved as the end
130+
#' of summer greendown date, and MidGreendown as the transition point of the
131+
#' first derivative of the autumn EVI2 curve.
129132
#' @param verbose logical. If `TRUE`, the progress will be reported.
130133
#'
131134
#' @return An object of class `BlspFit` will be returned. The object contains the
@@ -263,7 +266,8 @@ FitBLSP <- function(date_vec, vi_vec,
263266
p_li, mod,
264267
years, numYears,
265268
date_vec, vi_vec, weights_vec,
266-
cred_int_level
269+
cred_int_level,
270+
greendown_aware = opt$greendown_aware
267271
)
268272
} else {
269273
blsp_fit <- CalPhenoParam(

‎requirements/requirements.R‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,6 @@ pkgs <- c(
1313
"miniUI",
1414
"pkgdown",
1515
"devtools",
16-
"geojsonio",
17-
"geojsonR",
18-
"getPass",
1916
"minpack.lm",
2017
"rstac",
2118
"terra"

0 commit comments

Comments
 (0)