1. Introduction
In recent years, land surface temperature (LST) has become one of the most critical environmental covariates in scientific research [
1]. It plays an essential role in various fields such as urban remote sensing [
2,
3], vegetation analysis [
4], and cryosphere remote sensing [
5]. Using only conventional fixed-point observation on the ground makes obtaining the continuous spatial and temporal distribution of LST difficult, leaving satellite remote sensing the only feasible means of data acquisition [
6]. Remote sensing techniques, such as thermal infrared and microwave techniques, can be used to capture the LST. However, microwave sensors typically have lower spatial resolution for acquiring LST. The spatial and temporal resolution of thermal infrared LST detectors is mutually constrained, and weather often prevents complete data from being acquired. For instance, the MODIS sensor can acquire spatially continuous LST under cloud-free conditions, but the temporal continuity of LST is limited by the satellite revisit interval and cloud coverage.
Research has shown that over 60% of MODIS data in mid-latitude regions are af-fected by weather, with clouds being the main culprit [
7]. Compared with remote-sensing data, reanalysis data can offer spatiotemporally continuous surface temperature information. For example, the ERA5 reanalysis LST product from the European Centre for Medium-Range Weather Forecasts provides global hourly LST, albeit at a spatial resolution of only 0.1°. To satisfy the demand for LST with high spatiotemporal resolution, spatial downscaling of high temporal resolution LST has become an active research area in recent years [
8,
9]. The goal of such research is to improve the spatial resolution of downscaled LST.
Currently, LST downscaling methods can be broadly classified into two categories. The first category is based on image fusion, which involves fusing low spatial resolution and high temporal resolution LST data with high spatial resolution and low temporal resolution LST data to obtain high spatiotemporal resolution LST data [
10,
11]. Examples of such methods include the spatial and temporal adaptive reflectance fusion model (STARFM) [
12], the enhanced STARFM (ESTARFM) [
13], the spatiotemporal adaptive data fusion algorithm for temperature mapping [
10], the spatiotemporal integrated temperature fusion model [
11], and the deep learning based spatiotemporal temperature fusion network [
14]. These models all use MODIS LST data fusion technology and require only minimal auxiliary data while still maintaining good spatial texture for downscaled LST data. However, downscaled LST methods using image fusion techniques often overlook the radiative transfer of thermal infrared remote-sensing information [
15] and therefore lack clear physical mechanisms and scientific significance. Thus, to ensure the energy balance of surface thermal radiation information before and after downscaled LST, statistical regression downscaling methods based on the “scale-invariant relationship” hypothesis have emerged.
The second type of method for downscaling LST may be broadly classified as statistical regression based on scale factors. The principle of such methods is to establish a statistical regression relationship between LST and single or multiple land surface parameters and to assume that this statistical relationship is independent of spatial scales. Vegetation indices are the most used input parameters in statistical regression methods. For example, the DisTrad algorithm [
16] downscales LST by establishing a linear regression between LST and the normalized difference vegetation index (NDVI). The TsHARP algorithm [
17] downscales MODIS 1 km LST to 250 m and GOES 5 km LST to 1 km by establishing a linear relationship between LST and vegetation cover. The HUTS algorithm [
18] downscales LST by establishing a regression relationship between LST and the NDVI on the one hand and surface reflectance on the other hand. In addition, geographically weighted regression methods based on NDVI and digital elevation models have also been applied to downscaling research [
19]. Machine-learning methods can handle complex nonlinear statistical regression problems and have been applied with success to downscaling remote-sensing data. Previous studies have used artificial neural networks, random forests, support vector machines, and other methods combined with the NDVI, leaf area index, surface reflectance, digital elevation models, and other remote-sensing indices as input parameters to downscale LST [
20,
21,
22,
23]. Previous work evaluated the use of the machine-learning methods of random forests, neural networks, and support vector machines for downscaling different land-cover types and concluded that the random forest method is more stable and accurate than the other two methods [
24].
However, most input surface parameters (e.g., NDVI, leaf area index, FVC) for the above mentioned LST downscaling methods are calculated from visible or near-infrared remote-sensing images, which greatly limits the applicability of the downscaling methods under cloud cover. For all-weather acquisition of LST data, various schemes have been proposed for cloud-free LST reconstruction based on multi-source auxiliary information, LST inversion based on microwave remote-sensing data, and LST fusion based on synergistic thermal infrared and microwave remote-sensing data [
25,
26,
27]. However, whether these methods can be applied to the spatial downscaling of LST remains a critical issue that needs to be explored and addressed.
Given the discussion above, developing a method to downscale LST that is not limited by weather conditions and guarantees the spatial completeness and temporal continuity of the downscaled LST is an urgent problem that needs to be solved. The ERA5 reanalysis data from the European Centre for Medium-Range Weather Forecasts provides hourly LST with high temporal resolution, long time span, and easy accessibility. Previous studies show that the ERA5 LST product correlates strongly with actual LST, although the error is nonzero. However, whether it can be downscaled spatially to increase the spatiotemporal resolution requires further research. The present study thus proposes a pixel-wise temporal alignment iterative linear regression model (PTAILRM) for downscaling the ERA5 LST product based on the MODIS LST product and air temperature factors. The method involves time matching and iterative linear regression of each pixel, downscaling the continuously time-resolved, coarse-resolution ERA5 LST to a fine resolution and temporally sampling it to local time at whole-hour intervals. In this study, the ERA5 LST in the study area was downscaled to a resolution of 1000 m using the proposed method. The downscaling results of land-cover data and other representative methods were employed for qualitative evaluation, while the site-measured data were used for quantitative evaluation.
2. Study area and material
2.1 Study area
The study area encompasses Zhangye City and Haibei Tibetan Autonomous Prefecture, located in the middle reaches of the Heihe River Basin, with a total area of approximately 73 000 square kilometers. Zhangye City is in the middle section of the Hexi Corridor and is characterized by flat terrain dominated by cropland, grassland, bare soil, and wetland. The Haibei Tibetan Autonomous Prefecture is in the central part of the Qilian Mountains, with rugged terrain and elevations ranging from 2180 to 5287 meters above sea level. Areas above 3000 m account for more than 85% of the total prefectural area and are mainly covered by cropland, grassland, and bare soil.
Figure 1 shows the study area’s geographical location, elevation, and land cover.
2.2 Material
(i) The European Centre for Medium-Range Weather Forecasts provides the fifth-generation reanalysis climate dataset ERA5, which offers hourly reanalyzed data from 1979 onwards and is continuously updated. This dataset assimilates numerical weather forecasts, a significant amount of ground observation data, and satellite remote-sensing information and offers high temporal resolution and a long time span. Within the research area, we selected from ERA5-Land the hourly LST (LST) and air temperature (AT) products with a spatial resolution of 0.1° to conduct a long-term, all-weather analysis of LST downscaling.
Table 1 shows the data acquisition status.
(ii) The MOD11A1 and MYD11A1 products used in this study include daytime and nighttime surface temperature inversions from the Terra and Aqua satellites, respectively. The combined MOD11A1 and MYD11A1 products provide a sampling frequency of four times per day and are inverted by the generalized split-window algorithm to yield a surface temperature with approximately 1 km resolution [
28]. In addition, the study also uses the quality control (QC) layer and imaging time layer provided in the dataset.
Table 1 and
Table 2 give the data-acquisition status.
(iii) The observed meteorological data used to verify the accuracy were obtained from the National Tibetan Plateau Science Data Center.
Figure 2 shows the locations of the meteorological stations within the study area and
Table 3 provides information on the underlying surface type, altitude, and time range for each station. The observation dataset was the up-and-down longwave radiation recorded every 10 minutes by the meteorological stations [
29]. The observed surface temperature was calculated as follows by applying the Stefan–Boltzmann law to the up-and-down longwave radiation and MODIS emissivity dataset:
In Equation (1), L
↑ and L
↓ are the up and down longwave radiation, respectively, σ is the Stefan–Boltzmann constant [5.67×10−8 [W/(m
2·K
4)]]. The broadband emissivity ε
b is calculated according to Equation (2) by using the MODIS land surface emissivity dataset [
30], where ε
29, ε
30, and ε
31 are the narrowband emissivities of MODIS bands 29, 30, and 31, respectively.
3. Methodology
3.1 The PTAILRM downscaling method
3.1.1 Overview of PTAILRM
The PTAILRM downscaling method primarily consists of four components.
Step 1: Preprocessing of the ERA5 data and MODIS LST product, which involves spatial resampling of the ERA5 data and quality control of the MODIS LST product.
Step 2: Temporal alignment of the ERA5 data at each pixel, including conversion of each pixel’s Greenwich time to local time and interpolating the ERA5 product to match the transit time of the MODIS sensor.
Step 3: Pixel-wise iterative linear regression using the time series of ERA5 LST and AT as independent variables and the time series of MODIS LST at higher spatial resolution as the dependent variable. This step primarily involves establishing a linear regression model for each datum and iteratively determining the optimal time offset to obtain the regression coefficient matrix.
Step 4: Downsampling of the ERA5 LST and AT based on the coefficient matrix, followed by raster calculations to obtain downscaled LST. See
Figure 3 for a detailed description of the process.
3.1.2 Preprocessing of ERA5 and MODIS LST products
The ERA5 hourly raster data for LST and AT was resampled through bilinear interpolation to achieve 1000 m resolution that aligns with MODIS LST products. Although this process did not amplify the amount of information, empirical evidence shows that, compared with the nearest-neighbor interpolation method, this process mitigates the impact of significant gradient fluctuations at the edges of low-resolution pixels on the downscaling. The QC layer of the MOD11A1 and MYD11A1 LST products was used to mask regions affected by clouds, exclusively retaining pixels with average LST and emissivity errors below 1 K and 0.01, respectively. The aforementioned operation necessitates the use of the byte-type QC_Day and QC_Night data layers, specified in
Table 2, along with their corresponding bitmask information provided in
Table 4.
3.1.3 Pixel-wise temporal conversion and alignment
The ERA5 reanalysis data product is provided in Greenwich Mean Time, whereas the MODIS LST product is provided in local time. Due to the scanning characteristics of the MODIS instrument, variations exist in the local solar time of a given pixel on different revisit days [
31].
Figure 4 shows the local times corresponding to each MODIS overpass in the study area for 2012–2021, which are concentrated in four time periods (in decimal hours): 0.1–3, 9.9–11.9, 12.2–14.1, and 21.0–23.3.
LST correlates closely with the amount of solar radiation received on the ground, which in turn relates closely to the solar altitude, which varies with local time. Therefore, at the pixel scale, this study converts Greenwich Mean Time to local time, which is more closely related to changes in surface temperature than regional time. The time conversion formula is as follows:
where
,
, and
are the local time, Greenwich Mean Time, and longitude of pixel
i, respectively, with time measured in decimal hours.
After time conversion at the pixel level, we temporally align the ERA5 and MODIS products because the ERA5 hourly products are provided at hourly intervals, whereas the MODIS products have irregular sampling times. Aligning the two products temporally effectively minimizes errors caused by temporal inconsistencies, which is critical for maintaining the accuracy and reliability of the analysis.
If the function describing the LST as a function of time is continuous and differentiable everywhere, so interpolation methods can be used to obtain data at a given moment in a local time frame. The cubic spline function is commonly used for processing satellite time series data. Its primary advantages include continuous first and second derivatives, simplicity, stability, and high accuracy in calculations [
32,
33,
34,
35].
To address the problem of inconsistent sampling times between the ERA5 and MODIS LST, the high and continuous temporal resolution of the ERA5 product is used in combination with the view time data layer of the MODIS product to perform a cubic spline interpolation on a pixel-by-pixel basis. This temporally aligns the two products in the time domain (as illustrated in Figure 9).
3.1.4 Establishing the pixel-wise iterative linear regression model
As the periodic variation of solar radiation is the leading cause of variations in LST, the temporal patterns of LST for various land-cover types can be considered approximately synchronous. Given this, we hypothesize that the relationship between the coarse-resolution LST and fine-resolution LST for the same pixel after spatial and temporal registration is linear and remains stable within a specific time period, as expressed by the following formula:
In Equation (4), the subscript i signifies the pixel location, the independent variable t is the local time, ai and bi are the regression coefficients, Fi and Csi are the fine-resolution and coarse-resolution LSTs, respectively, and is the residual. Notably, a time offset parameter tsi is introduced in the equation to mitigate errors arising from asynchronous temperature changes between Fi and Ci. It is evident that tsi is relatively small and fluctuates around zero, and its value varies depending on the location of pixel i.
Given that the value selected for
tsi should not be large, an iterative regression can be done with a given step size within a time window centered at zero to determine the
tsi that corresponds to the regression with the highest R
2. Assume the relationship between
tsi and R
2 is expressed by
where ω is the window size, τ
si is the value of
tsi in the time window, and R
2si is the R
2 corresponding to τ
si. By taking specific steps to vary
tsi within the given time window and performing iterative regression,
tsi can be determined by finding the maximum R
2si, as given by Equation (6). Furthermore, land and atmospheric interactions, along with energy exchange, significantly affect LST and result in a strong correlation between LST and AT [
36,
37,
38]. To further improve the precision of the downscaling process, we introduce an AT factor and establish the following linear relationship between coarse-resolution AT and fine-resolution LST:
Equation (7) is like Equation (4) and uses the subscript i to indicate the pixel location, along with an independent variable t and the regression coefficients aai and bai. Additionally, Equation (7) includes fine-resolution LST (Fi) and coarse-resolution AT (Cai) as well as residual and time offset tai, which can be acquired in the same way as done for tsi.
In accordance with Equations (4) and (7), the fine-resolution LST is
and the simulated fine-resolution LST is
In Equations (8) and (9), the regression coefficients are ai, bi, and ci, the residual is , and i is the simulated fine-resolution LST. The other variables have already been introduced.
Equation (9) is a relatively simple model comprising only three unknown parameters: ai, bi, and ci, which can be determined via linear regression. tsi and tai can be obtained through iterative regression, with a time window ω of one hour used herein. In practical downscaling operations, t is expressed in cumulative decimal hours, a unit also employed in the view time layer of MODIS LST products, and its values are determined by the cloud-free imaging time series of MODIS LST products within a specific time period, which is one month in this study, ensuring sufficient data for regression analysis.
Upon iteratively performing linear regression to determine the parameter matrix, it can be combined with the temporally and spatially aligned time series of ERA5 LST and AT for band operations, ultimately yielding the hourly LST at a spatial resolution of 1000 m. The temporal coverage extends from January 1, 2012 at 00:00:00 to December 31, 2021 at 23:00:00.
3.2 Validation strategy
3.2.1 Strategy for qualitative validation
The strategy for qualitative validation primarily encompasses two aspects. First, it entails evaluating the spatiotemporal rationality of LST by presenting maps of it for representative periods and analyzing daily and seasonal LST variations in conjunction with land-cover types. Second, it involves comparisons with other representative downscaling methods.
We compare the ESTARFM [
13] with the proposed PTAILRM. ESTARFM functions by initially inputting two sets of coarse-resolution and fine-resolution images from the same or neighboring time periods, followed by inputting a coarse-resolution image acquired at the target date, ultimately yielding a fine-resolution image at the target date. The merit of this fusion technique resides in its capacity to account for the spatial proximity between adjacent and target pixels, spectral discrepancies, and temporal fluctuations while augmenting the accuracy of fusion outcomes in regions with pronounced heterogeneity.
Identifying two pairs of suitable ERA5 and MODIS LSTs is relatively easy and allows us to downscale low-resolution ERA5 LST acquired at any given time within a proximate time range, thereby facilitating comparison with the downscaled results PTAILRM.
3.2.2 Quantitative validation strategy
This investigation precisely evaluates the MODIS overpass instances under cloud-free conditions by using goodness-of-fit indices, namely, R² and the mean absolute error (MAE). This approach verifies only the accuracy of the downscaled LST during imaging periods when the sky was transparent, yet it can be amalgamated with land-cover classifications, digital elevation models, and cloud-free pixel count (CPC) to scrutinize how disparate surface factors affect the accuracy. The emphasis is placed on examining how spatial factors affect the accuracy.
The methodology for assessing accuracy based on data from meteorological stations is used extensively with remote-sensing LST products [
39,
40,
41,
42]. This approach is constrained by the quantity of data from meteorological stations and the relatively homogeneous land-cover types surrounding these stations; however, it most important feature is the temporal resolution, which facilitates the comparison and analysis of downscaled accuracy over an entire time period without being limited by weather conditions. To avoid the influence of scale factors, data from meteorological stations must serve as a singular data source for both input and accuracy validation [
43]. This method underscores the relationship between temporal factors and downscaling error. This study thus adheres to the following principles:
(1) Data from meteorological stations are used as a single input source for analyzing the accuracy of downscaling models. That is, ERA5 surface temperature is downscaled separately by using the PTAILRM at the spatial scale of the infrared observatory of the five meteorological stations, which prevents scale effects from influencing the accuracy analysis.
(2) The time sampling used for the regression input data is consistent with the MODIS surface-temperature product’s clear-sky transit time on the pixel scale corresponding to each meteorological station. In other words, the sampling strategy simulates the MODIS transit time.
4. Downscaling results
4.1 Qualitative validation of downscaling results
4.1.1 Validation of the spatiotemporal regularity of downscaled LST
Figure 6 shows the downscaled LST on the first day of representative months for each season in 2021. To accommodate spatial constraints, a 4-hour temporal interval was adopted. The LST derived from the PTAILRM downscaling has rich texture, pronounced heterogeneity, superior spatial and temporal integrity, and few missing pixels.
This study uses a monthly average method to calculate and represent the hourly diurnal variations of the downscaled LST for three distinctive surface-cover types: urban regions, vegetation, and ice- and snow-covered landscapes (see
Figure 7). Correspondingly, their deviations from daily mean temperatures and daily maximum temperatures are tabulated to facilitate a comprehensive analysis and evaluation (see
Figure 8).
As shown in
Figure 7, the diurnal variation of the downscaled LST reveals a gradual rise in LST within the study area from sunrise (approximately at 6:00) to its zenith around 14:00 local time. Subsequently, between 14:00 and 18:00, the LST gradually declines. Throughout the nighttime period (from 18:00 to 6:00 the following day), the temperature remains relatively low with predictable fluctuations.
The overall trend shown in
Figure 8 suggests that the daily average temperature and the maximum temperature difference decrease in the following order: urban built-up areas, vegetation, and snow and ice-covered areas; With regards to the downscaled LST seasonal variation pattern, the highest daily average LST occurs during summer (June–August), whereas the lowest average LST occurs during winter (December–February). During spring (March–May), the average temperature is slightly higher than during autumn (September–November).
4.1.2 Comparison of PTAILRM with Alternative Downscaling Model
Figure 9 shows the original ERA5 LST at 2 pm on the first day of representative months for each season and the LST downscaled by using PTAILRM and ESTARFM. The reason we selected 2 pm is that the LST typically reaches its daily maximum around this time and exhibits maximal spatial texture.
Figure 9 shows that both downscaling models increase the spatial texture details of the LST vis-à-vis the original LST. However, the ESTARFM downscaled LST has more missing pixels than the PTAILRM downscaled LST. This phenomenon occurs at all four selected time points, especially on January 1 and October 1. Additionally, the downscaling of ESTARFM is not stable. For January 1, its overall temperature is lower, and the spatial distribution of LST deviates significantly from the original LST.
Moreover, temperature anomalies occur at the lower part of the image. For January 1, the PTAILRM downscaled LST better preserves the spatial distribution of the original LST. For April 1, ESTARFM’s downscaled image is noisy and contains serious spatial discontinuities, where the PTAILRM avoids these issues. For July 1 and October 1, the ESTARFM downscaled results are better, with fewer missing pixels. However, the spatial texture details of the ESTARFM downscaled LST differ significantly from those of the PTAILRM in areas with high and low LSTs.
4.2 Quantitative validation of downscaled results
4.2.1 Validation of MODIS transit times under cloud-free conditions
This section analyzes the accuracy of the PTAILRM downscaling results at the monthly scale in combination with factors such as CPC, elevation, and land-cover type.
Figure 10 shows the spatial distribution of the coefficient of determination (R
2) and MAE produced by the PTAILRM and CPC. This model produces a high goodness-of-fit for each month, with mean R
2 and MAE values reaching 0.87 and 2.7 K, respectively. However, the results for April and June are less precise than the other months, with R
2 = 0.77 and 0.83 and MAE = 3.5 and 3.4 K, respectively. Furthermore, the mean CPC values of 55 and 53 are lower during these time periods. The difference in precision between the remaining months is not significant, with mean R
2 and CPC values exceeding 0.83 and 64 and MAE values below 3.0 K.
Figure 11(a)–11(c) show the statistical distribution of MAE and CPC for each month, as well as the relationship between the distribution of MAE data and the variation of CPC. As the CPC increases, the fitting accuracy remains stable, whereas the uncertainty increases as the amount of data decreases.
Figure 11(d) and 11(e) show the statistical fluctuation of MAE with respect to elevation and land-cover type. The shape and pattern of the box plot indicate that the MAE varies less in terms of elevation and land-cover type when compared to CPC. In
Figure 11(d), the MAE distribution has high values in the middle elevations (around 4000 m) and low values at extreme elevations. In
Figure 11(d), the overall MAE for cropland and urban areas is relatively low, whereas, for grassland and wetland, it is relatively high compared with other land-cover types.
4.2.2 Validation for all time periods based on data from meteorological stations
Figure 12 shows the statistical distribution of errors in the downscaled results on an hourly scale for annual mean, using error, mean error (ME), and MAE as evaluation metrics.
Figure 12(a
1)–12(a
5) and 12(b
1)–12(b
5) correspond to clear and non-clear sky periods for the five meteorological stations, respectively, while
Figure 12(a6) and 12(b6) show the mean statistical distribution of various indicators for the five meteorological stations under cloud-free and overcast conditions, respectively.
Figure 13 shows the statistical analysis of the ME and MAE for each time period, along with the absolute differences between each pair of time periods.
As shown in
Figure 12, the ME and MAE values under cloud-free and overcast conditions exhibit greater fluctuations and relatively higher values during the daytime (from approximately 6 am to 6 pm). In contrast, during the nighttime, the fluctuations are smaller, and the values are lower. Additionally,
Figure 12(a
6) and 12(b
6) show the varied ME under cloud-free and overcast conditions, respectively. Generally, ME values during cloud-free conditions are negative in the morning and positive in the afternoon, implying that downscaled LST tends to be overestimated in the morning and underestimated in the afternoon. Under overcast conditions, ME values are negative in the morning, indicating overestimation, and fluctuate near zero in other time periods.
Figure 13 indicates that, during MODIS transit periods, cloud-free periods, and nighttime periods, the absolute values of both the ME and MAE are less than those during non-transit periods, overcast periods, and daytime periods. Notably, the ME for cloud-free and overcast conditions indicate an overall underestimation and overestimation, respectively, of the downscaled LST. The minimum and maximum MEs occur during the overcast and cloud-free time periods, with values of −0.98 and 0.13 K, respectively. Conversely, the maximum and minimum MAEs occur during daytime and nighttime segments, with values of 2.01 and 0.85 K, respectively. The largest absolute difference in MEs occurs between cloud-free and overcast segments, with a difference of 1.11 K, whereas in MAEs it occurs between daytime and nighttime segments, with a difference of 1.16 K.
5. Discussion
5.1 Qualitative performance of PTAILRM
Using the PTAILRM to downscale the LST produces high spatiotemporal completeness, rich spatial texture, and compliance with various spatiotemporal variations of surface cover. These results demonstrate that the proposed method is suitable for downscaling the ERA5 LST product with high temporal resolution and coarse spatial resolution. Additionally, the pixel-wise temporal registration and iterative regression ensure high spatial heterogeneity of the downscaled results.
Compared with the ESTARFM, the PTAILRM results are more stable and complete, which is attributed to the assumption of linear variation in land surface reflectance over time made by ESTARFM and to the impact of missing pixels and the requirement on input data of high spatial completeness. The proposed method, however, does not require input data with high spatial completeness.
The analysis of spatial factors contributing to errors shows that CPC is more sensitive to MAE than surface-cover types and altitude because the quantity of input data affects the stability of iterative regression. Unstable fitting occurs for grassland and wetland, which may be attributed to frequent changes in grassland conditions due to local animal husbandry and changes in soil moisture content in wetlands. The MAE is relatively low for extreme altitudes (i.e., high and low) but relatively high for moderate altitudes, possibly due to the frequent changes in snowline in the moderate-altitude Qilian Mountains in the study area [
44], which would decrease in the accuracy of the model. In contrast, surface-cover changes beyond the snowline range are relatively stable, resulting in more accurate fitting.
5.2 Quantitative accuracy of PTAILRM
The analysis based on data from meteorological stations reveals that the errors during non-view times, overcast conditions, and daytime are greater than those during view times, cloud-free conditions, and nighttime. View and cloud-free times produce more accurate results than non-view and overcast times because the downscaled LST lacks samples from the latter. The error during nighttime is lower than that during daytime because the temperature changes more dramatically during the day, which increases the gap between the simulated LST and the measured LST. Additionally, under clear-sky conditions, the downscaled LST is overestimated in the morning and underestimated in the afternoon, which may be attributed to the morphological differences between the downscaled LST time series and the original LST time series. Under overcast conditions, the downscaled LST tends to be overestimated, possibly because the data used for iterative regression all came from remote-sensing products acquired under clear-sky conditions.
The maximum absolute difference in the ME occurs between cloud-free and overcast conditions, indicating that the ME is sensitive to these two conditions because the downscaled LST is generally underestimated and overestimated under cloud-free and overcast conditions, respectively. The maximum absolute difference in the MAE occurs between daytime and nighttime, indicating that the MAE is more responsive to changes in solar radiation than to other conditions, such as when the satellite passes or if the weather is clear. In contrast, this phenomenon does not occur in the pairwise comparison of other time periods.
6. Conclusions
The PTAILRM is suitable for downscaling the long-term all-weather ERA5 product. It uses only clear-sky MODIS LST as input data, and the method of iterative regression by pixel-wise temporal registration ensures a relatively high accuracy of downscaled LST, even for cloudy periods. Additionally, the PTAILRM produces accurate LST maps and is well-suited for spatiotemporal analysis of LST for various land-cover types. Furthermore, given that this study considered a region characterized by complex terrain and harsh weather conditions, it is reasonable to expect that greater downscaling accuracy could be achieved in regions with gentle topography and mild climates.
PTAILRM assumes a pixel-wise linear correlation between temporal LST at fine resolution and temporal LST and AT at coarse spatial resolution. This correlation includes an unknown time offset variable, which can be determined using iterative linear regression. The validity of these assumptions is confirmed by quantitatively evaluating the accuracy of the downscaled LST. Moreover, the model is relatively simple with few fitting parameters, so it requires only a small amount of remote-sensing LST data as input, making it easy to apply. This simple model avoids the overfitting that is common of complex models and which can cause excessive image noise.
However, the method proposed in this study still has the following limitations:
The method requires long-term coarse-resolution and fine-resolution data as input to ensure the stability of iterative regression. In other words, lower CPC values may lead to instability in the iterative regression process, which could result in larger errors in the downscaling of LST.
Due to the mutual constraints of the temporal and spatial resolution of remote-sensing surface temperature products, only MODIS data were selected as the data source in this study because they have relatively high temporal and spatial resolution. Other remote-sensing surface-temperature products with higher spatial resolution usually have lower temporal resolution, resulting in uneven time sampling and poor downscaled results. These limitations lead to a lack of comparative analysis of the downscaled accuracy of this method for different resolutions.
The iterative regression method in the PTAILRM partly eliminates errors inherent in the ERA5 LST product. However, how these errors, along with the downscaling model errors, affect the overall error requires further research.
Author Contributions
N.W. proposed the theory of the PTAILRM and designed the experiment. N.W. and J.T. performed the overall data analysis and drafted the manuscript. S.S. and Q.T. provided comments and suggestions for the manuscript and checked the writing. All authors have read and agreed to the published version of the manuscript.
Funding
This research was funded by Open Fund of State Key Laboratory of Urban and Re-gional Ecology (grant number: SKLURE2023-2-6), State Key Laboratory of Remote Sensing Science (grant number: OFSLRSS202119), National Natural Science Foundation of China (grant number: 42101321), and China Postdoctoral Science Foundation (grant number: 2021M701653).
Data Availability Statement
Acknowledgments
We express our gratitude to the following organizations for graciously providing the data utilized in this study: the European Centre for Medium-Range Weather Forecasts (ECMWF), for providing ERA5 surface temperature, the Land Processes Distributed Active Archive Center (LPDAAC) for providing MOD11A1 and MYD11A1 products, and the National Tibetan Plateau Science Data Center (TPBC) for providing in-situ measurements data. Moreover, the authors wish to thank the anonymous reviewers for their valuable comments, which considerably contributed to the quality of the paper, and to Prof. Wenfeng Zhan for valuable discussions and advice.
Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper.
References
- Wu, J.; Xia, L.; On Chan, T.; Awange, J.; Zhong, B. Downscaling Land Surface Temperature: A Framework Based on Geographically and Temporally Neural Network Weighted Autoregressive Model with Spatio-Temporal Fused Scaling Factors. ISPRS J. Photogramm. Remote Sens. 2022, 187, 259–272. [Google Scholar] [CrossRef]
- Hrisko, J.; Ramamurthy, P.; Yu, Y.; Yu, P.; Melecio-Vázquez, D. Urban Air Temperature Model Using GOES-16 LST and a Diurnal Regressive Neural Network Algorithm. Remote Sens. Environ. 2020, 237, 111495. [Google Scholar] [CrossRef]
- Wang, D.; Chen, Y.; Hu, L.; Voogt, J.A.; Gastellu-Etchegorry, J.-P.; Krayenhoff, E.S. Modeling the Angular Effect of MODIS LST in Urban Areas: A Case Study of Toulouse, France. Remote Sens. Environ. 2021, 257, 112361. [Google Scholar] [CrossRef]
- Julien, Y.; Sobrino, J.A. The Yearly Land Cover Dynamics (YLCD) Method: An Analysis of Global Vegetation from NDVI and LST Parameters. Remote Sens. Environ. 2009, 113, 329–334. [Google Scholar] [CrossRef]
- Langer, M.; Westermann, S.; Boike, J. Spatial and Temporal Variations of Summer Surface Temperatures of Wet Polygonal Tundra in Siberia - Implications for MODIS LST Based Permafrost Monitoring. Remote Sens. Environ. 2010, 114, 2059–2069. [Google Scholar] [CrossRef]
- Dash, P.; Göttsche, F.-M.; Olesen, F.-S.; Fischer, H. Land Surface Temperature and Emissivity Estimation from Passive Sensor Data: Theory and Practice-Current Trends. Int. J. Remote Sens. 2002, 23, 2563–2594. [Google Scholar] [CrossRef]
- Shui-sen Chen; Xiu-zhi Chen; Wei-qi Chen; Yong-xian Su; Da Li A Simple Retrieval Method of Land Surface Temperature from AMSR-E Passive Microwave Data—A Case Study over Southern China during the Strong Snow Disaster of 2008. Int. J. Appl. Earth Obs. Geoinf. 2011, 13, 140–151. [CrossRef]
- Weng, Q.; Fu, P. Modeling Diurnal Land Temperature Cycles over Los Angeles Using Downscaled GOES Imagery. ISPRS J. Photogramm. Remote Sens. 2014, 97, 78–88. [Google Scholar] [CrossRef]
- Zhan, W.; Chen, Y.; Zhou, J.; Wang, J.; Liu, W.; Voogt, J.; Zhu, X.; Quan, J.; Li, J. Disaggregation of Remotely Sensed Land Surface Temperature: Literature Survey, Taxonomy, Issues, and Caveats. Remote Sens. Environ. 2013, 131, 119–139. [Google Scholar] [CrossRef]
- Weng, Q.; Fu, P.; Gao, F. Generating Daily Land Surface Temperature at Landsat Resolution by Fusing Landsat and MODIS Data. Remote Sens. Environ. 2014, 145, 55–67. [Google Scholar] [CrossRef]
- Wu, P.; Shen, H.; Zhang, L.; Goettsche, F.-M. Integrated Fusion of Multi-Scale Polar-Orbiting and Geostationary Satellite Observations for the Mapping of High Spatial and Temporal Resolution Land Surface Temperature. Remote Sens. Environ. 2015, 156, 169–181. [Google Scholar] [CrossRef]
- Feng Gao; Masek, J.; Schwaller, M.; Hall, F. On the Blending of the Landsat and MODIS Surface Reflectance: Predicting Daily Landsat Surface Reflectance. IEEE Trans. Geosci. Remote Sensing 2006, 44, 2207–2218. [CrossRef]
- Zhu, X.; Chen, J.; Gao, F.; Chen, X.; Masek, J.G. An Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model for Complex Heterogeneous Regions. Remote Sens. Environ. 2010, 114, 2610–2623. [Google Scholar] [CrossRef]
- Yin, Z.; Wu, P.; Foody, G.M.; Wu, Y.; Liu, Z.; Du, Y.; Ling, F. Spatiotemporal Fusion of Land Surface Temperature Based on a Convolutional Neural Network. IEEE Trans. Geosci. Remote Sens. 2021, 59, 1808–1822. [Google Scholar] [CrossRef]
- Zhan, W.; Chen, Y.; Zhou, J.; Li, J.; Liu, W. Sharpening Thermal Imageries: A Generalized Theoretical Framework From an Assimilation Perspective. IEEE Trans. Geosci. Remote Sens. 2011, 49, 773–789. [Google Scholar] [CrossRef]
- Kustas, W.; Norman, J.; Anderson, M.; French, A. Estimating Subpixel Surface Temperatures and Energy Fluxes from the Vegetation Index-Radiometric Temperature Relationship. Remote Sens. Environ. 2003, 85, 429–440. [Google Scholar] [CrossRef]
- Agam, N.; Kustas, W.P.; Anderson, M.C.; Li, F.; Neale, C.M.U. A Vegetation Index Based Technique for Spatial Sharpening of Thermal Imagery. Remote Sens. Environ. 2007, 107, 545–558. [Google Scholar] [CrossRef]
- Dominguez, A.; Kleissl, J.; Luvall, J.C.; Rickman, D.L. High-Resolution Urban Thermal Sharpener (HUTS). Remote Sens. Environ. 2011, 115, 1772–1780. [Google Scholar] [CrossRef]
- Duan, S.-B.; Li, Z.-L. Spatial Downscaling of MODIS Land Surface Temperatures Using Geographically Weighted Regression: Case Study in Northern China. IEEE Trans. Geosci. Remote Sens. 2016, 54, 6458–6469. [Google Scholar] [CrossRef]
- Hutengs, C.; Vohland, M. Downscaling Land Surface Temperatures at Regional Scales with Random Forest Regression. Remote Sens. Environ. 2016, 178, 127–141. [Google Scholar] [CrossRef]
- Keramitsoglou, I.; Kiranoudis, C.T.; Weng, Q. Downscaling Geostationary Land Surface Temperature Imagery for Urban Analysis. IEEE Geosci. Remote Sens. Lett. 2013, 10, 1253–1257. [Google Scholar] [CrossRef]
- Sismanidis, P.; Keramitsoglou, I.; Kiranoudis, C.T.; Bechtel, B. Assessing the Capability of a Downscaled Urban Land Surface Temperature Time Series to Reproduce the Spatiotemporal Features of the Original Data. Remote Sens. 2016, 8. [Google Scholar] [CrossRef]
- Yang, G.; Pu, R.; Huang, W.; Wang, J.; Zhao, C. A Novel Method to Estimate Subpixel Temperature by Fusing Solar-Reflective and Thermal-Infrared Remote-Sensing Data With an Artificial Neural Network. IEEE Trans. Geosci. Remote Sens. 2010, 48, 2170–2178. [Google Scholar] [CrossRef]
- Li, W.; Ni, L.; Li, Z.-L.; Duan, S.-B.; Wu, H. Evaluation of Machine Learning Algorithms in Spatial Downscaling of MODIS Land Surface Temperature. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2019, 12, 2299–2307. [Google Scholar] [CrossRef]
- Duan, S.-B.; Li, Z.-L.; Leng, P. A Framework for the Retrieval of All-Weather Land Surface Temperature at a High Spatial Resolution from Polar-Orbiting Thermal Infrared and Passive Microwave Data. Remote Sens. Environ. 2017, 195, 107–117. [Google Scholar] [CrossRef]
- Huang, C.; Duan, S.-B.; Jiang, X.-G.; Han, X.-J.; Leng, P.; Gao, M.-F.; Li, Z.-L. A Physically Based Algorithm for Retrieving Land Surface Temperature under Cloudy Conditions from AMSR2 Passive Microwave Measurements. Int. J. Remote Sens. 2019, 40, 1828–1843. [Google Scholar] [CrossRef]
- Zhao, G.; Gao, H.; Cai, X. Estimating Lake Temperature Profile and Evaporation Losses by Leveraging MODIS LST Data. Remote Sens. Environ. 2020, 251, 112104. [Google Scholar] [CrossRef]
- Wan, Z.; Dozier, J. A Generalized Split-Window Algorithm for Retrieving Land-Surface Temperature from Space. IEEE Trans. Geosci. Remote Sens. 1996, 34, 892–905. [Google Scholar] [CrossRef]
- Liu, S.; Li, X.; Xu, Z.; Che, T.; Xiao, Q.; Ma, M.; Liu, Q.; Jin, R.; Guo, J.; Wang, L.; et al. The Heihe Integrated Observatory Network: A Basin-Scale Land Surface Processes Observatory in China. Vadose Zone J. 2018, 17, 180072. [Google Scholar] [CrossRef]
- Wang, K.; Wan, Z.; Wang, P.; Sparrow, M.; Liu, J.; Zhou, X.; Haginoya, S. Estimation of Surface Long Wave Radiation and Broadband Emissivity Using Moderate Resolution Imaging Spectroradiometer (MODIS) Land Surface Temperature/Emissivity Products. J. Geophys. Res.: Atmos. 2005, 110. [Google Scholar] [CrossRef]
- Sharifnezhadazizi, Z.; Norouzi, H.; Prakash, S.; Beale, C.; Khanbilvardi, R. A Global Analysis of Land Surface Temperature Diurnal Cycle Using MODIS Observations. J. Appl. Meteorol. Climatol. 2019, 58, 1279–1291. [Google Scholar] [CrossRef]
- Mao, F.; Li, X.; Du, H.; Zhou, G.; Han, N.; Xu, X.; Liu, Y.; Chen, L.; Cui, L. Comparison of Two Data Assimilation Methods for Improving MODIS LAI Time Series for Bamboo Forests. Remote Sens. 2017, 9, 401. [Google Scholar] [CrossRef]
- Singh, R.K.; Liu, S.; Tieszen, L.L.; Suyker, A.E.; Verma, S.B. Estimating Seasonal Evapotranspiration from Temporal Satellite Images. Irrig. Sci. 2012, 30, 303–313. [Google Scholar] [CrossRef]
- Wüst, S.; Wendt, V.; Linz, R.; Bittner, M. Smoothing Data Series by Means of Cubic Splines: Quality of Approximation and Introduction of a Repeating Spline Approach. Atmos. Meas. Tech. 2017, 10, 3453–3462. [Google Scholar] [CrossRef]
- Zhang, H.; Pu, R.; Liu, X. A New Image Processing Procedure Integrating PCI-RPC and ArcGIS-Spline Tools to Improve the Orthorectification Accuracy of High-Resolution Satellite Imagery. Remote Sens. 2016, 8, 827. [Google Scholar] [CrossRef]
- Benali, A.; Carvalho, A.C.; Nunes, J.P.; Carvalhais, N.; Santos, A. Estimating Air Surface Temperature in Portugal Using MODIS LST Data. Remote Sens. Environ. 2012, 124, 108–121. [Google Scholar] [CrossRef]
- Mostovoy, G.V.; King, R.L.; Reddy, K.R.; Kakani, V.G.; Filippova, M.G. Statistical Estimation of Daily Maximum and Minimum Air Temperatures from MODIS LST Data over the State of Mississippi. Gisci. Remote Sens. 2006, 43, 78–110. [Google Scholar] [CrossRef]
- Yang, Y.Z.; Cai, W.H.; Yang, J. Evaluation of MODIS Land Surface Temperature Data to Estimate Near-Surface Air Temperature in Northeast China. Remote Sens. 2017, 9. [Google Scholar] [CrossRef]
- Ermida, S.L.; Soares, P.; Mantas, V.; Göttsche, F.-M.; Trigo, I.F. Google Earth Engine Open-Source Code for Land Surface Temperature Estimation from the Landsat Series. Remote Sens. 2020, 12, 1471. [Google Scholar] [CrossRef]
- Ermida, S.L.; Trigo, I.F.; DaCamara, C.C.; Göttsche, F.M.; Olesen, F.S.; Hulley, G. Validation of Remotely Sensed Surface Temperature over an Oak Woodland Landscape — The Problem of Viewing and Illumination Geometries. Remote Sens. Environ. 2014, 148, 16–27. [Google Scholar] [CrossRef]
- Li, H.; Li, R.; Yang, Y.; Cao, B.; Bian, Z.; Hu, T.; Du, Y.; Sun, L.; Liu, Q. Temperature-Based and Radiance-Based Validation of the Collection 6 MYD11 and MYD21 Land Surface Temperature Products Over Barren Surfaces in Northwestern China. IEEE Trans. Geosci. Remote Sens. 2021, 59, 1794–1807. [Google Scholar] [CrossRef]
- Yu, W.; Ma, M.; Yang, H.; Tan, J.; Li, X. Supplement of the Radiance-Based Method to Validate Satellite-Derived Land Surface Temperature Products over Heterogeneous Land Surfaces. Remote Sens. Environ. 2019, 230, 111188. [Google Scholar] [CrossRef]
- Hong, F.; Zhan, W.; Göttsche, F.-M.; Lai, J.; Liu, Z.; Hu, L.; Fu, P.; Huang, F.; Li, J.; Li, H.; et al. A Simple yet Robust Framework to Estimate Accurate Daily Mean Land Surface Temperature from Thermal Observations of Tandem Polar Orbiters. Remote Sens. Environ. 2021, 264, 112612. [Google Scholar] [CrossRef]
- Guo, Z.; Wang, N.; Shen, B.; Gu, Z.; Wu, Y.; Chen, A. Recent Spatiotemporal Trends in Glacier Snowline Altitude at the End of the Melt Season in the Qilian Mountains, China. Remote Sens. 2021, 13, 4935. [Google Scholar] [CrossRef]
Figure 1.
Geographical location and land-cover types in the study area. Panel (a) shows the geography and elevation of the study area, and panel (b) shows the primary land cover types.
Figure 1.
Geographical location and land-cover types in the study area. Panel (a) shows the geography and elevation of the study area, and panel (b) shows the primary land cover types.
Figure 2.
Location of meteorological stations within the study area. The map displays a standard pseudo color image of the study area; detailed information on stations 1–5 is available in Table 3.
Figure 2.
Location of meteorological stations within the study area. The map displays a standard pseudo color image of the study area; detailed information on stations 1–5 is available in Table 3.
Figure 3.
Flow chart of PTAILRM. The green rounded rectangles in the figure represent input data or data generated during the processing, the blue containers show data operations, and the plus sign indicates the combined input of data.
Figure 3.
Flow chart of PTAILRM. The green rounded rectangles in the figure represent input data or data generated during the processing, the blue containers show data operations, and the plus sign indicates the combined input of data.
Figure 4.
MODIS view time for 2012–2021.
Figure 4.
MODIS view time for 2012–2021.
Figure 5.
Spline interpolation of time series data.
Figure 5.
Spline interpolation of time series data.
Figure 6.
Spatial representation of downscaled LST. Shown are the downscaled LST for the first day of January, April, July, and October to represent each season in 2021. The results are displayed at 4-hour intervals.
Figure 6.
Spatial representation of downscaled LST. Shown are the downscaled LST for the first day of January, April, July, and October to represent each season in 2021. The results are displayed at 4-hour intervals.
Figure 7.
Monthly average of hourly LST between 2012–2021 for representative land covers.
Figure 7.
Monthly average of hourly LST between 2012–2021 for representative land covers.
Figure 8.
Bar graphs of daily mean temperature and maximum temperature difference between 2012–2021 for three major land-cover types.
Figure 8.
Bar graphs of daily mean temperature and maximum temperature difference between 2012–2021 for three major land-cover types.
Figure 9.
Comparison of downscaled LST produced by PTAILRM and ESTARFM. (a)–(d) Downscaled LSTs at 2:00 pm on January 1, April 1, July 1, and October 1 of 2021, respectively. The dotted ellipses show areas of significant difference between the two methods. The black regions indicate missing data at the corresponding pixels.
Figure 9.
Comparison of downscaled LST produced by PTAILRM and ESTARFM. (a)–(d) Downscaled LSTs at 2:00 pm on January 1, April 1, July 1, and October 1 of 2021, respectively. The dotted ellipses show areas of significant difference between the two methods. The black regions indicate missing data at the corresponding pixels.
Figure 10.
Spatial distribution of R2, MAE, and CPC for each month between 2012–2021.
Figure 10.
Spatial distribution of R2, MAE, and CPC for each month between 2012–2021.
Figure 11.
Pixel-wise MAE distribution during cloud-free MODIS views at 1000 m scale from 2012 to 2021.
Figure 11.
Pixel-wise MAE distribution during cloud-free MODIS views at 1000 m scale from 2012 to 2021.
Figure 12.
Error, ME, and MAE statistical distribution on hourly scale from 2012 to 2021. The gray areas represent the sampling periods for the input data used in the regression, which correspond to the overpass times of MODIS.
Figure 12.
Error, ME, and MAE statistical distribution on hourly scale from 2012 to 2021. The gray areas represent the sampling periods for the input data used in the regression, which correspond to the overpass times of MODIS.
Figure 13.
ME and MAE for each time period from 2012 to 2021. The absolute differences between the ME and MAE values for view time versus non-view time, cloud-free time versus overcast time, and daytime versus nighttime are 1, 2, and 3, respectively.
Figure 13.
ME and MAE for each time period from 2012 to 2021. The absolute differences between the ME and MAE values for view time versus non-view time, cloud-free time versus overcast time, and daytime versus nighttime are 1, 2, and 3, respectively.
Table 1.
Characteristics of data used in this study.
Table 1.
Characteristics of data used in this study.
Dataset |
Temporal range |
Spatial resolution |
Temporal resolution |
ERA5 LST |
2012/01/01–2021/12/31 |
0.1° |
1 h |
ERA5 AT |
2012/01/01–2021/12/31 |
0.1° |
1 h |
MOD11A1 |
2012/01/01–2021/12/31 |
1000 m |
12 h |
MYD11A1 |
2012/01/01–2021/12/31 |
1000 m |
12 h |
Table 2.
MODIS LST product used in this study.
Table 2.
MODIS LST product used in this study.
Data Layer |
Units |
Description |
LST_Day_1km |
Kelvin |
Daytime LST |
LST_Night_1km |
Kelvin |
Nighttime LST |
Day_view_time |
Hours |
Local time of day observation |
Night_view_time |
Hours |
Local time of night observation |
QC_Day |
- |
Daytime LST Quality Indicators |
QC_Night |
- |
Nighttime LST Quality indicators |
Table 3.
Information for each meteorological station.
Table 3.
Information for each meteorological station.
ID |
Site name |
Altitude |
Sampling time interval |
Time range |
1 |
Zhangye |
1460 m |
10 min |
2012/01/01–2021/12/31 |
2 |
Yakou |
4148 m |
10 min |
2012/01/01–2021/12/31 |
3 |
Jingyangling |
3750 m |
10 min |
2012/01/01–2021/12/31 |
4 |
Huazhaizi |
1731 m |
10 min |
2012/01/01–2021/12/31 |
5 |
Dashalong |
3739 m |
10 min |
2012/01/01–2021/12/31 |
Table 4.
Bitmask description for QC_Day and QC_Night layer.
Table 4.
Bitmask description for QC_Day and QC_Night layer.
Bits location |
Value |
Description |
Flag name |
Bits 0 & 1 |
0 |
LST produced, good quality, not necessary to examine more detailed QA |
Mandatory QA flags |
1 |
LST produced, other quality, recommend examination of more detailed QA |
2 |
LST not produced due to cloud effects |
3 |
LST not produced primarily due to reasons other than cloud |
Bits 2 & 3 |
0 |
Good data quality |
Data quality flag |
1 |
Other quality data |
2 |
TBD |
3 |
TBD |
Bits 4 & 5 |
0 |
Average emissivity error ≤ 0.01 |
Emissivity error flag |
1 |
Average emissivity error ≤ 0.02 |
2 |
Average emissivity error ≤ 0.04 |
3 |
Average emissivity error > 0.04 |
Bits 6 & 7 |
0 |
Average LST error ≤ 1K |
LST error flag |
1 |
Average LST error ≤ 2K |
2 |
Average LST error ≤ 3K |
3 |
Average LST error > 3K |
|
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).