Downscaling MODIS land surface temperature to 90 m using random forest regression to assess transferability

MODIS spectral bands and LST

MOD09GA.061 dataset (https://doi.org/10.5067/MODIS/MOD09GA.061) was used to select six spectral bands, namely blue, green, red, near infrared (NIR), shortwave infrared 1 (SWIR1), and shortwave infrared 2 (SWIR2). MOD11A1.061 dataset (https://doi.org/10.5067/MODIS/MOD11A1.061) was used to derive the LST. MODIS LST products store raw data as scaled integers to optimize storage efficiency and maintain precision. A scaling factor of 0.02 is applied to convert these integer values into physically meaningful LST values in Kelvin (K). Finally, temperatures in degrees Celsius (°C) are derived by subtracting 273.15 from the Kelvin results.

MODIS vegetation indices

The MOD09GA.061 dataset was used to derive the following vegetation indices:

Normalized difference vegetation index

(1) $$NDVI={\displaystyle \frac{NIR-Red}{NIR+Red}}.$$

Normalized difference water index

(2) $$NDWI={\displaystyle \frac{Green-NIR}{Green+NIR}.}$$

Modified normalized difference water index

(3) $$MNDWI={\displaystyle \frac{Green-SWIR1}{Green+SWIR2}.}$$

Normalized difference built-up index

(4) $$NDBI={\displaystyle \frac{SWIR1-NIR}{SWIR1+NIR}.}$$

Soil-adjusted vegetation index

(5) $$SAVI={\displaystyle \frac{\left(NIR-Red\right)\left(1+L\right)}{NIR+Red+L}.}$$

Here, $L$ is the soil brightness correction factor.

Vegetation indices are integral to LST downscaling as they provide essential information about vegetation cover, influencing surface temperatures (Wang et al., 2021; Pandey et al., 2023).

For digital elevation model data, two scene IDs, P5_PAN_CS_N22_000_E085_000_30m and P5_PAN_CS_N23_000_E085_000_30m from the PAN sensor of the CartoSat-1 satellite (https://bhoonidhi.nrsc.gov.in/bhoonidhi/index.html) were mosaiced to encompass the regions of interest.

Landsat 9 LST retrieval

The NDVI threshold-based emissivity method is one of the approaches to retrieve Landsat 9 LST.

The raw sensor-specific quantized and calibrated digital numbers (DNs) stored in Landsat imagery (https://earthexplorer.usgs.gov/) are converted into physically quantifiable radiance units.

(6) $$L_{\lambda }=\left(M_L\times Q_{cal}\right)+A_L.$$Here, $M_L$ and $A_L$ are the band-specific multiplicative (gain) and additive (offset) rescaling factors, respectively, provided in the imagery’s metadata file, and $Q_{cal}$ denotes the quantized and calibrated pixel values (DNs).

The inverse of Planck’s function is applied to convert radiance into brightness temperature.

(7) $$T_b={\displaystyle \frac{K_2}{\ln \left({\displaystyle \frac{K_1}{L_{\lambda }}+1}\right)}.}$$Here, $K_1$ and $K_2$ are band-specific thermal conversion constants extracted from the metadata, and $L_{\lambda }$ is the spectral radiance in $W{m}^{-2}s{r}^{-1}\mu m^{-1}$.

To adjust the brightness temperature for actual surfaces, each pixel’s surface emittance of thermal radiation is considered. The emissivity depends on surface composition, as vegetation, soil, rock, water, and urban infrastructure emit radiation differently. To quantify this, $NDVI$ is computed. The following equation calculates the proportion of vegetation, $P_v$ within a pixel, representing the extent of green vegetation cover to quantify surface heterogeneity, using NDVI extremes for “complete vegetation” ( $NDV{I}_v$) and “bare soil” ( $NDV{I}_s$).

(8) $$P_v={\left({\displaystyle \frac{NDVI-NDV{I}_v}{NDV{I}_v-NDV{I}_s}}\right)}^2.$$Higher vegetation cover has higher emissivity, while bare soil or rock has lower emissivity. Using $P_v$ and known emissivity values for dense vegetation and bare soil, land surface emissivity (LSE) for each pixel is dynamically calculated using Eq. (9).

(9) $$\epsilon =\{\begin{array}{cc}\hfill {\epsilon }_{s\lambda },& NDVI<NDV{I}_s\hfill \\ {\epsilon }_{v\lambda }{P}_v+{\epsilon }_{s\lambda }\left(1-P_v\right)+C,& NDV{I}_s\le NDVI\le NDV{I}_v\hfill \\ \hfill {\epsilon }_{v\lambda },& NDVI>NDV{I}_v\hfill \end{array}.$$Here, ${\mathrm{\epsilon }}_{s\lambda }$ and ${\mathrm{\epsilon }}_{v\lambda }$ are the emissivities of bare soil and pure vegetation, respectively, and $C$ represents surface roughness, taken as a constant.

Empirical NDVI-based emissivity relationships have been used extensively in the literature. These methods ensure that LSE is not treated as a constant but is spatially dynamic and tailored to the surface conditions (Sobrino, Jiménez-Muñoz & Paolini, 2004; Momeni & Saradjian, 2007; Sobrino et al., 2008; Li et al., 2013). For a smooth transition, emissivity was modeled as ε = 0.004P_v + 0.986 in this study.

With LSE determined, the at-sensor brightness temperature is converted to LST by accounting for the reduced recorded thermal radiance caused by the surface’s deviation from an ideal blackbody behavior.

(10) $$LST=\frac{T_b}{1+\left(\frac{\lambda T_b}{\rho }\right)\ln \left(\epsilon \right)}.$$Here, $T_b$ is the brightness temperature, $\lambda$ is the center wavelength of the TIR band (band 10) in $\mu m$, and $\rho$ is $1.4388\times {10}^{-2}mK$ derived from Eq. (11), and $\mathrm{\epsilon }$ is the LSE.

(11) $$\rho ={\displaystyle \frac{hc}{\sigma }}.$$Here, $h$ is the Planck’s constant ( $6.626\times {10}^{-34}Js$), σ is the Boltzmann constant ( $1.38\times {10}^{-23}J{K}^{-1}$), and $c$ is the velocity of light ( $2.998\times {10}^{8}m{s}^{-1}$).

Model construction using RF

The RF model is a powerful ensemble machine-learning method introduced by Breiman (2001). It is a low-bias, low-variance model as it tends to generalize well and perform consistently across training and test datasets. RF effectively handles large datasets that exhibit complex relationships. It is capable of noise reduction and outlier removal, and avoids overfitting problems by averaging the outputs of many decorrelated trees (Prasad et al., 2020; Njuki, Mannaerts & Su, 2020). The decorrelation is achieved through bootstrapping and selecting a random subset of features at each split. The interpretability through feature importance is a significant advantage.

There are different feature selection methods, namely Boruta, RFE, and multicollinearity methods, which have been applied in various studies (Prasad et al., 2023; Abdo et al., 2024). This study initially considered thirteen independent features for the primary region. To mitigate the curse of dimensionality and prevent overfitting, an RFECV method using gradient boosting was employed, selecting eight relevant independent features to enhance the model’s generalization capability (Barzani et al., 2024). The eight features include four spectral bands (blue, near infrared (NIR), shortwave infrared (SWIR)1, and SWIR2), two vegetation indices (MNDWI and SAVI), MODIS LST, and digital elevation model (DEM) processed from two CartoSat-1 scenes as the only topographic feature. Finer resolution (90 m) Landsat 9-retrieved LST was used as the target variable. The 90-m scale was fixed to preserve the band’s native resolution maximally. The Landsat LST data was processed and directly exported to a 90-m scale in Google Earth Engine (GEE). All the MODIS variables and DEM were resampled to 90 m using bilinear interpolation for spatial congruence with the target variable.

The model (Fig. 2) was built using optimal parameters obtained after randomly exploring 50 combinations from the grid of values of five parameters (Table 3), namely n_estimators, criterion, max_depth, min_samples_leaf, and max_features, with a 5-fold cross-validation using the RandomizedSearchCV method of the open-source scikit-learn Python library to minimize the bias-variance trade-off for improving the model’s predictive accuracy and achieve robustness (Bergstra & Bengio, 2012). The minimum value in the set of values for min_samples_leaf was carefully chosen to be 10 to restrain the growth of the decision trees to the deepest level, avoiding the risk of fitting any potential outliers. The parameter n_jobs was set to −1 to involve all the CPU cores for parallelization. The RandomizedSearchCV method was employed to reduce computational time while exploring diverse combinations. All the computing was done by an AMD Ryzen 5 5600H processor (3.30 GHz) with 16 GB RAM (3,200 MHz) in a Windows 11 environment.

To evaluate the model’s generalization capability and avoid overfitting, the dataset was partitioned into training and test subsets in the ratio 7:3, ensuring the model was validated on unseen data before applying it to the secondary region. Once trained, the RF model was used to predict the LST for the secondary region, enhancing the spatial resolution of the MODIS LST. The following metrics (Eqs. (12)–(14)) were used to evaluate the downscaling performance of the RF regression model for both the primary and secondary regions. Since the model was trained on 70% of the data from the primary region, the remaining 30% was used to evaluate the model to obtain an unbiased estimate of its predictive performance. The trained model was then applied to the secondary region for the reproducibility of similar results observed in the primary region’s test set.

(12) $$\begin{array}{c}MAE={\displaystyle \frac{1}{n}\sum _{i=1}^n\left|LS{T}_{a,i}-LS{T}_{p,i}\right|}\end{array}$$

(13) $$\begin{array}{c}RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{\left(LS{T}_{a,i}-LS{T}_{p,i}\right)}^2}}\end{array}.$$Here, $LS{T}_{a,i}$ and $LS{T}_{p,i}$ are actual and predicted LST values for $i$-th pixel, respectively, and $n$ is the number of pixels.

(14) $$\begin{array}{c}{R}^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{\left(LS{T}_{a,i}-LS{T}_{p,i}\right)}^2}{{\sum }_{i=1}^n{\left(LS{T}_{a,i}-\overline{LS{T}_a}\right)}^2}}\end{array}.$$Here, $\overline{LS{T}_a}$ is the mean around the actual LST values, and $n$ is the number of pixels.

Feature importance by RF

Determining the importance of features is essential in LST, natural resources, and hazard mapping (Prasad et al., 2022, 2024, 2025). The feature importance plots of the primary region (Fig. 3) weigh heavily on three factors, namely, SWIR1, MODIS LST, and DEM. MODIS LST contributed 68.89% and 63.04% in January and February, respectively, to the model’s predictions. In these months, SWIR1 was the next highest contributor, influencing the predictions by 13.46% and 16.93%, respectively. The influence of MODIS LST in May and June elevated to 84.77% and 93.10%, respectively, undermining all other features. A contrasting dynamic was observed in December, where the importance of SWIR1 dominated all other features, influencing 58.02% of the model’s predictions. There are notable spikes in the importance of the other factors, namely NIR, blue, SWIR2, SAVI, MNDWI, and DEM, in November compared to January, February, May, and June. The most important feature for spatial downscaling was MODIS LST, followed by SWIR1, DEM, SWIR2, and MNDWI. In contrast, the blue band, SAVI, and NIR band were among the less important factors.

Downscaling performance

The predictions on the test data across all the months aligned closely with the actual LST values (Table 4). The mean temperatures were consistent, with trifling deviations in the months of February and December. The model performed exceptionally well in May and June, as evident by the root mean square error (RMSE) and R² values. In June, the RMSE was the lowest, and the R² value maximally justified the 1:1 ideal line, explaining 96% of the variability. The spread of LST values, reflecting the level of information acquired (Mukherjee, Joshi & Garg, 2017; Maithani et al., 2022), was significantly closer to the actual variance in June. The model predicted a variance of 6.83 against the actual 7.25. This month’s minimum and maximum predicted LST values of 27.10 °C and 40.30 °C depicted the closest resemblance to the actual LST range, justifying the actual variance. The predictions were weaker in the months of February and December, with lower variances. The R² was relatively lower (0.76) in November and December, indicating the model’s inefficiency in fully capturing the variability. The test predictions also revealed a general tendency to compress the maximum and marginally elevate the minimum temperatures, affecting the variability. Substantial compression in the temperature ranges can be observed in the months of February and November, with deviations of 3.56 °C and 3.50 °C. The model’s tendency to underestimate high temperatures is offset by its simultaneous overestimation of low temperatures, keeping the mean temperatures consistent with the actual values. This predictive behavior is a fundamental characteristic of RF regression and is attributed to the model’s approach to averaging the outputs from all the decision trees.

The scatterplots revealed a substantial concentration of predictions against the actual LST values on the 1:1 ideal line (Fig. 4). The model underestimated the higher temperatures in all the months except May and June. There was no considerable disarray of predictions in May and June, and a small number of predictions were affected by a high error magnitude. Overall, the magnitude of errors was relatively greater on the edges of the actual temperature spectrum across all the months.

The error histograms exhibited the model’s unbiased predictability, as the errors were symmetrically distributed around zero for all the months (Fig. 5). The bell-shaped curves suggest a normal distribution of errors, with most errors accumulating within the range of −2 to +2 °C. A high percentage of errors has a ±1 °C offset from the ideal value of 0 across all the months. The narrower distributions in May and June, particularly June, highlighted the model’s strong performance with fewer errors and more accurate predictions.

Model transferability

The trained RF model was applied to the secondary region to assess its transferability. It showed a moderate drop in efficiency, reflecting the challenges of applying models to unknown spatial contexts. The secondary region analyses yielded R² values between 0.65 and 0.76, with RMSE values ranging from 0.85 °C to 1.18 °C (Table 5). The predictions in January had an RMSE of 1.01 and an R² of 0.76. The variability was underrepresented in the predicted temperatures, with the predicted variance (2.87) lower than the actual variance (4.22). A similar pattern was observed in February with an RMSE of 0.96 and an R² of 0.74. The predicted mean LST (22.34 °C) closely matched the actual mean (22.36 °C). Contrary to the test results, the model’s performance declined considerably in May and June in the secondary region with RMSEs of 1.18 and 1.01, respectively. The lower R² values indicated a weaker alignment of predictions against the actual values. With an RMSE of 0.85 and an R² value of 0.74, the November predictions captured trends to a large extent and maintained variability. In December, the model yielded an RMSE of 1.12 and an R² of 0.73. The predicted variance (3.10) in this month was significantly lower than the actual variance (4.74), reflecting limitations in capturing local variability. The model’s performance in simulating the temperature ranges varied seasonally within the secondary region. While the model depicted fair agreement with the observed temperature ranges during the peak warm months (especially the month of May), it consistently underestimated the temperature variations in the transitional months (January, February, November, and December), mostly estimating narrower ranges than observed.

The secondary region scatterplots (Fig. 6) showed a derangement of some predictions, suggesting weaker correlations in all the months, notably the May and June months. There is a consistent underestimation of higher temperature ranges in all the months. While most predictions are tightly clustered along the 1:1 ideal line, indicating smaller errors, larger errors tend to cluster at the extremes of the target range, suggesting the model’s sensitivity in predicting extreme values. Since the model was trained in the primary region, the clusters of incorrectly predicted points in these plots indicate that the secondary region exhibits different spatial features that were never represented in the training data. Despite variations in R² values, the model maintained consistent performance throughout the six temporal points, never dropping below 0.65, suggesting fundamental reliability.

All six error histograms (Fig. 7) show a single, sharply peaked (high kurtosis) distribution around zero, especially in February, November, and December, indicating that most predictions align with the actual values. A slight offset in all the plots, evident from the right elongated tails in the distributions, further validates the underestimation in the predictions. The tails beyond ±2 °C for all the months confirm the existence of larger errors, but the rapid drop-off in frequency indicates they are relatively rare compared to the strong concentrations near zero.

The secondary region’s actual and predicted 90-m LST maps are presented in Fig. 8. The temperature spectrum is represented in a color gradient from blue to red, with blue indicating lower temperatures and red indicating higher temperatures. A visual comparative analysis of the maps depicts considerable similarity, making the RF model suitable for downscaling.

The considerable white sparse areas in the predicted per-pixel RMSE maps (Fig. 9) depict close alignment with the actual LST values with lesser deviation (≤1 °C). Most predictions can be observed within a derangement of 2 °C, with a significantly low number of predicted temperature values exceeding an RMSE beyond 3 °C.

Comparative evaluation of RF against XGBoost and CatBoost

To contextualize the performance of the RF model, a comparative analysis was done using XGBoost and categorical boosting (CatBoost). These models were trained and evaluated using the exact temporal points, feature sets, datasets, and validation procedures to ensure an appropriate comparison of their downscaling capabilities.

The XGBoost model, introduced by Chen & Guestrin (2016), mitigates the common overfitting issue in traditional boosting models by incorporating regularization techniques. These penalties help control model complexity, preventing overfitting. The model’s objective function is optimized using a second-order Taylor expansion for faster convergence.

In the context of regression, CatBoost, introduced by Prokhorenkova et al. (2018) is a powerful gradient boosting algorithm that predicts a continuous target variable by minimizing the default RMSE loss function. CatBoost trains each tree on a sequence of randomly shuffled datasets, ensuring that the error for a particular data point is evaluated on a model that has not been trained on that specific data point. This leads to more robust models that generalize better to unseen data.

To ensure optimal performance, hyperparameter tuning was employed for both algorithms using a grid search strategy similar to the RF-based methodology. The number of different hyperparameter combinations (n_iter) was set to 50, keeping it consistent with the RF model. Tables 6 and 7 specify the relevant parameters that were tuned with the grid of values and the optimal combinations.

The performance statistics for both primary and secondary regions are specified in Tables 8 and 9. The XGBoost model produced a performance profile highly similar to the RF model on the primary region’s test set (Table 8). Across all the months, the metrics for both models were almost identical, suggesting that XGBoost is an equally feasible model for this dataset, offering a similar predictive advantage. However, in the secondary region, it consistently underperformed compared to RF (Table 9). XGBoost registered a higher RMSE and a lower R² value than the RF model across all six temporal points (Fig. 10). This suggests that the XGBoost model was slightly more biased to the primary region’s data and therefore less generalizable than RF.

As illustrated in Fig. 11, the residual distribution for the XGBoost model in the secondary region is broadly comparable to the RF model. However, a key distinction appears in the November and December plots, where the RF model exhibits elongated right tails. These tails signify large-magnitude errors that are notably absent in the XGBoost results. Despite this, the RF model achieved a lower overall RMSE because of its higher volume of extremely accurate (near-zero error) predictions, effectively outweighing the statistical penalty from the infrequent outliers.

In contrast, the CatBoost model consistently straggled behind the RF model across all the evaluated months and metrics on the primary region’s test set (Table 8). CatBoost produced higher mean absolute error (MAE) and RMSE values and lower R² scores than RF in all the months. In May, CatBoost’s RMSE and R² were 0.75 and 0.88 respectively, both notably worse than the RF’s RMSE of 0.65 and R² of 0.91. This pattern of weak performance was consistent, as observed in December, where CatBoost’s RMSE and R² were 0.79 and 0.71 respectively, compared to the RF’s more accurate RMSE (0.72) and R² (0.76). In the secondary region, the CatBoost model managed to match or even slightly outperform RF in a few months, e.g., in February, when its RMSE was 0.95 (Table 9) compared to RF’s 0.96 (Table 5) and November, when its RMSE was 0.84 compared to RF’s 0.85. But it experienced a pronounced failure in June. In this month, its R² declined to just 0.32 (Fig. 12), whereas the RF model maintained an R² of 0.65 (Fig. 6). This extreme inconsistency makes CatBoost unstable and therefore, inapposite for downscaling as a generalized model based on this study.

The CatBoost residuals plots (Fig. 13) exhibit close similarity with XGBoost in the months of January and May. Wider distributions can be observed in CatBoost in the months of June and December compared to XGBoost. Though inconsiderable in frequency, residuals of magnitude greater than 5 °C and extending up to 10 °C can be observed in the month of November in CatBoost. All the CatBoost residual plots show wider distributions for all the months, contrasting with the RF plots that depict strong concentration of residuals within ±2 °C.