Downscaling land surface temperature data using edge detection and block-diagonal Gaussian process regression

Accurate and high-resolution estimation of land surface temperature (LST) is crucial in estimating evapotranspiration, a measure of plant water use and a central quantity in agricultural applications. In this work, we develop a novel statistical method for downscaling LST data obtained from NASA’s ECOSTRESS mission, using high-resolution data from the Landsat 8 mission as a proxy for modeling agricultural field structure. Using the Landsat data, we identify the boundaries of agricultural fields through edge detection techniques, allowing us to capture the inherent block structure present in the spatial domain. We propose a block-diagonal Gaussian process (BDGP) model that captures the spatial structure of the agricultural fields, leverages independence of LST across fields for computational tractability, and accounts for the change of support present in ECOSTRESS observations. We use the resulting BDGP model to perform Gaussian process regression and obtain high-resolution estimates of LST from ECOSTRESS data, along with uncertainty quantification. Our results demonstrate the practicality of the proposed method in producing reliable high-resolution LST estimates, with potential applications in agriculture, urban planning, and climate studies.

💡 Research Summary

The paper presents a novel statistical framework for downscaling land surface temperature (LST) observations from NASA’s ECOSTRESS mission (≈70 m resolution) to a higher spatial resolution (30 m) by exploiting the fine‑scale structural information contained in Landsat 8 imagery. The authors first generate an accurate partition of the study area into agricultural fields using edge detection on an averaged Landsat image. They evaluate several edge detectors and find that Meta’s Segment Anything Model (SAM) combined with a custom post‑processing step yields the most reliable, non‑overlapping field boundaries, which define a set of blocks (R_1,\dots,R_N).

Within each block the mean LST temporal behavior is modeled with a simple Fourier regression: an annual sinusoid for Landsat (because of its sun‑synchronous orbit) and both annual and diurnal sinusoids for ECOSTRESS (which provides multiple observations per day). The Fourier coefficients are estimated by ordinary least squares, and the residuals are assumed to follow a zero‑mean Gaussian process (GP).

To capture spatial dependence while preserving computational tractability, the authors introduce a block‑diagonal Gaussian process (BDGP). For block (r) the covariance is the standard squared‑exponential kernel with its own variance (\sigma_r^2) and length‑scale (\ell_r). The overall covariance matrix is block‑diagonal, implying independence between fields. This assumption is justified by the physical expectation that temperature fields are only weakly correlated across field boundaries and dramatically reduces the computational cost from (\mathcal{O}(n^3)) to a sum of (\mathcal{O}(n_r^3)) operations, where (n_r) is the number of pixels in block (r).

The observation model explicitly accounts for the “change of support” problem: ECOSTRESS measurements are modeled as the true LST field convolved with a Gaussian blur kernel (derived from the instrument’s point‑spread function, with (\sigma_{\text{blur}}\approx0.97) pixels) plus additive sensor noise ((\sigma_{\text{sensor}}\approx0.1) K). Landsat observations are treated as unblurred. To keep the model identifiable, the authors assume the spatial length‑scale is the same for both sensors ((\ell_{\text{ES}}=\ell_{\text{LS}})).

Parameter estimation proceeds in two stages. With the block partition fixed, the authors fit (\sigma_{\text{LS}}) and (\ell_{\text{LS}}) for each block by maximizing the likelihood of the Landsat residuals, using four temporally well‑separated dates to reduce temporal correlation. A small nugget term ensures numerical stability. Next, using the same block structure, they estimate (\sigma_{\text{ES}}) for ECOSTRESS by a method‑of‑moments formula that subtracts sensor noise and accounts for the blur kernel; the final estimate is the median across the same four dates.

For the actual downscaling, the fitted BDGP for ECOSTRESS is used in a kriging step on a 30 m grid. Because blurring introduces weak dependence across block borders, the authors enlarge each block by a 4‑σ blur neighbourhood ((R’r)) and perform kriging independently within these extended regions. The kriging predictor is the standard linear combination ( \hat f{\text{ES}}(x^*) = k^{\top}(K_{R’r}+\sigma{\text{sensor}}^2 I)^{-1} y_{R’_r}), and the associated prediction variance is computed analytically, providing a full uncertainty map.

The methodology is demonstrated over the Salton Sea region in Southern California. The SAM‑derived field boundaries align well with visible agricultural patterns. Fourier fits capture the dominant seasonal cycle, and the residuals exhibit spatial structure that is effectively modeled by the BDGP. Compared to recent machine‑learning downscaling approaches (random forests, deep neural networks), the proposed method reduces oversmoothing at field edges and delivers calibrated uncertainty estimates. Computationally, the block‑wise formulation enables parallel processing and makes the approach feasible for large remote‑sensing datasets.

In conclusion, the paper contributes (1) an automated, high‑quality edge‑detection pipeline for defining spatial blocks, (2) a block‑diagonal GP that balances model fidelity with scalability, (3) an explicit treatment of sensor blur and noise within a Bayesian framework, and (4) a complete downscaling solution that delivers both high‑resolution LST maps and pixel‑wise uncertainty. The authors suggest that the framework can be extended to other remotely sensed variables that exhibit block‑like spatial heterogeneity (e.g., soil moisture, air quality) and that future work will address dynamic block structures, incorporation of additional ancillary data, and real‑time operational deployment.