A Bayesian model to estimate land surface phenology parameters with harmonized Landsat 8 and Sentinel-2 images

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. Three case studies are presented and used to explore aspects of the proposed model. The first case study, conducted over forested pixels within an 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 case study, conducted over the same pixels using only 2018 data, compares annual LSP parameter estimates from our proposed models with those generated using methods described in Bolton et al. (2020). The third case study, conducted on a small area of interest within the HLS tile on an annual time-step (2014–2019), further examines the impact of sample size and choice of likelihood on annual LSP parameter estimates in addition to assessing potential for the proposed models to inform LSP change detection analysis. 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. The 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 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. Modeling functions in the rsBayes package and supplementary code sections are threaded, allowing for the use of multiple processing cores to further speed up model fitting for massive datasets. Using a 64 CPU workstation, the first case study analysis took ~3 days to run using the Beta likelihood model. However, processing time decreases linearly as the number of CPU cores increases. We expect that run times for this LSP modeling approach will decrease substantially as the power of new computing systems increases over time.