The Landsat satellite data applications have been broadened by the improvement in data consistency, e.g., advances in calibration and atmospheric correction algorithms. However, there are still inter-sensor inconsistencies of surface reflectance caused by differences in sensor spectral response and residuals in calibration and atmospheric correction. Such inconsistencies tend to entangle with the land surface change, especially for the subtle changes such as vegetation responses to climate change. We propose a Bayesian approach to reduce inter-sensor inconsistencies by making Landsat 5, 7 and 8 reflectance consistent with MODIS nadir BRDF-adjusted reflectance (NBAR). The Bayesian MODIS-consistent Landsat reflectance should be similar to the MODIS NBAR when it is degraded to the MODIS resolution considering the point spread function (PSF). The MODIS NBAR PSF is modelled as a Gaussian function with empirically derived and band specific parameters. The relationships between the MODIS-consistent Landsat reflectance and the two available reflectance (i.e., the MODIS NBAR and observed Landsat reflectance) are built in the Bayesian framework. The Bayesian MODIS-consistent Landsat reflectance is then derived using the computationally efficient conjugate gradient descent method. The algorithm was applied to Landsat 5, 7 and 8 surface reflectance over five study areas in Brazil, China, France, Russia and USA, and compared with a linear regression-based reflectance adjustment method. Visual difference was evident in the MODIS NBAR and the observed Landsat reflectance and the relative root-mean-squared differences (rRMSD) ranged from 3.16% to 67.72% and R² from 0.037 to 0.934. Both the Bayesian and linearly MODIS-consistent Landsat reflectance are visually more consistent with the MODIS NBAR. The average rRMSD and R² improvements of the Bayesian MODIS-consistent Landsat reflectance over the observed Landsat reflectance were 21.36% and 0.19, 3.38% and 0.11, 5.25% and 0.07, 3.14% and 0.11, 3.49% and 0.08, and 5.72% and 0.09 for the blue, green, red, near infrared, and two shortwave infrared bands, respectively. The linearly MODIS-consistent Landsat reflectance had less rRMSD improvements: 18.26%, 1.46%, 3.02%, 1.86%, 1.97% and 3.36% for the six bands respectively, and had no R² improvements. The linear regression derived using the whole image pixels has low performance in adjusting some extreme reflectance values that are not well represented in the image. In contrast, the Bayesian method derives the reflectance image that is consistent to MODIS at pixel level.