  We develop a Bayesian Land Surface Phenology (LSP) model and examine its
performance using Enhanced Vegetation Index (EVI) observations derived from the
Harmonized Landsat Sentinel-2 (HLS) dataset. Building on previous work, we
propose a double logistic function that, once couched within a Bayesian model,
yields posterior distributions for all LSP parameters. We assess the efficacy
of the Normal, Truncated Normal, and Beta likelihoods to deliver robust LSP
parameter estimates. Two case studies are presented and used to explore aspects
of the proposed model. The first, conducted over forested pixels within a HLS
tile, explores choice of likelihood and space-time varying HLS data
availability for long-term average LSP parameter point and uncertainty
estimation. The second, conducted on a small area of interest within the HLS
tile on an annual time-step, further examines the impact of sample size and
choice of likelihood on LSP parameter estimates. Results indicate that while
the Truncated Normal and Beta likelihoods are theoretically preferable when the
vegetation index is bounded, all three likelihoods performed similarly when the
number of index observations is sufficiently large and values are not near the
index bounds. Both case studies demonstrate how pixel-level LSP parameter
posterior distributions can be used to propagate uncertainty through subsequent
analysis. As a companion to this article, we provide an open-source \R package
\pkg{rsBayes} and supplementary data and code used to reproduce the analysis
results. The proposed model specification and software implementation delivers
computationally efficient, statistically robust, and inferentially rich LSP
parameter posterior distributions at the pixel-level across massive raster time
series datasets.
