Land surface phenology (LSP) is a consistent and sensitive indicator of climate change effects on Earth’s vegetation. Existing methods of estimating LSP require time series densities that, until recently, have only been available from coarse spatial resolution imagery such as MODIS (500 m) and AVHRR (1 km). LSP products from these datasets have improved our understanding of phenological change at the global scale, especially over the MODIS era (2001-present). Nevertheless, these products may obscure important finer scale spatial patterns and longer-term changes. Therefore, we have developed a Bayesian hierarchical model to retrieve complete annual sequences of LSP from Landsat imagery (1984-present), which has medium spatial resolution (30 m) but relatively sparse temporal frequency. Our approach uses Markov Chain Monte Carlo (MCMC) sampling to quantify individual phenometric uncertainty, which is especially important when considering long time series with variable observation quality and density, but has rarely been demonstrated. The estimated spring LSP had strong agreement with ground phenology records at Harvard Forest (R2 = 0.87) and Hubbard Brook Experimental Forest (R2 = 0.67). The estimated LSP were consistent with the recently released 30 m LSP product, MSLSP30NA, in its time period of 2016 to 2018 (R2 of 0.86 and 0.73 for spring and autumn phenology, respectively). Our Bayesian hierarchical model is an important step forward in extending medium resolution LSP records back in time as it accomplishes both critical goals of retrieving annual LSP from sparse time series and accurately estimating uncertainty.