Extracting Land Surface Albedo from Landsat 9 Data in GEE Platform to Support Climate Change Analysis

: Land surface albedo is a relevant variable in many climatic, environmental, and hydrological studies; its monitoring allows researchers to identify changes on the Earth’s surface. The open satellite data that is provided by the USGS/NASA Landsat mission is quite suitable for estimating this parameter through the remote sensing technique. The purpose of this paper is to evaluate the potentialities of the new Landsat 9 data for retrieving Earth’s albedo by applying da Silva et al.’s algorithm (developed in 2016 for the Landsat 8 data) using the Google Earth Engine cloud platform and R software. Two urban areas in Southern Italy with similar geomorphologic and climatic characteristics were chosen as study sites. After obtaining thematic maps of the albedos here, a statistical analysis and comparison among the Landsat 8 and Landsat 9 results was performed considering the entire study areas and each land use/land cover class that is provided by the Copernicus Urban Atlas 2018. This approach was also applied to the data after being filtered through Tukey’s test (used to detect and remove outliers). The analysis showed a very good correlation between the Landsat 8 and Landsat 9 estimations ( ρ > 0.94 for both sites), with some exceptions that were related to some mis-corresponding values. Furthermore, the Landsat 8 and Landsat 9 outliers were generally overlapping. In conclusion, da Silva et al.’s approach appears to also be reasonably applicable to the Landsat 9 data despite some radiometric differences.


Introduction
Land surface albedo (usually called "albedo") is a key geophysical variable that has attracted the interest of many researchers who are studying climate change, Earth's surface energy budget, and water resource management.Indeed, albedo information (which quantifies a surface's ability to reflect incident solar radiation) allows for an evapotranspiration (ET) estimate, which is particularly valuable for water resource management and tackling the issue of drought [1][2][3][4][5].In urban areas, lower ET values correspond to impermeable surfaces, implying that the heat from the Sun is stored and later released into the atmosphere, resulting in an urban heat island (UHI) effect [5][6][7][8][9][10][11][12][13][14][15].Monitoring surface albedo is, therefore, relevant for evaluating the impacts of climate changes, which are defined in the United Nations' (UN) 17 Sustainable Development Goals (SDGs) for safeguarding the Earth from global threats [3,13,[16][17][18].
Open Earth Observation (EO) satellite data plays an important role in predicting environmental variables and tracking changes on the Earth's surface at different scales of investigation.While low-resolution data is mainly used for global monitoring, medium-resolution Landsat and Sentinel-2 images (with geometric resolutions of 30 and 10 m, respectively) are adequate for regional and city scale research [1,4,10,16,[19][20][21].Because assessing the consequences of climate change necessitates the data spanning of a long time period, the multidecadal data that is provided by the Landsat mission (which has been in operation since 1972) is usually chosen over Sentinel-2 data (which has only been available since 2015) [10][11][12][13].
The Landsat satellites that are currently in orbit include Landsat 8 (launched in 2013) and Landsat 9 (in space since 2021).The sensors of Landsat 8's Operational Land Imager (OLI) and Landsat 9's OLI-2 have comparable characteristics, such as geometric, temporal, and spectral resolutions.This responds to the need to keep satellites monitoring natural resources and processes, thereby identifying long-term changes on the Earth's surface.Nonetheless, the Landsat 9 instrument has a higher radiometric resolution (14 bits) than the OLI sensor (which has a 12-bit resolution).As a result, the debut of the Landsat 9 platform provides still-not-much-explored opportunities for environmental monitoring from EO data [2,10,[22][23][24][25][26][27].
There are various approaches mentioned in the literature for determining albedos from remote sensing (RS) images.Bidirectional reflectance distribution function (BRDF) angular models and the narrow-to-broadband conversion methods are the most common.The former models are more complex and require more parameters than the latter ones, which are commonly adopted to simplify the estimation procedure [1,3].In relation to the latter, one of the most-often-used is that of Liang [28], which is an easy-to-use technique that is suited for a wide range of EO data.In 2008, Tasumi et al. [29] suggested a more complicated strategy for managing Landsat 5 and 7 data based on atmospheric correction and a radiative transfer model.Da Silva et al. (2016) [30] introduced the first algorithm for Landsat 8 OLI images, which enabled the retrieval of broadband albedos from top-of-atmosphere (TOA) images.This approach has been proven to produce accurate estimates despite its need for a significant amount of implementation and operating time (owing to the necessity to calculate atmospheric transmissivity in order to apply the atmospheric adjustment) [10,16,30].
To speed up RS image acquisition and processing times, a novel free cloud-computing-based platform (Google Earth Engine [GEE]) has been recently released to efficiently store, process, and analyze geospatial big data.In fact, this engine is made up of a high-performance computational service that is linked to an interactive programming interface as well as a multi-petabyte data catalog that includes a massive number of raw and preprocessed open satellite data sets.Because of the benefits that were discussed above, GEE is a powerful tool with more functionalities than standard geographic information system (GIS) software [10,[31][32][33][34][35][36][37].However, it is frequently paired with the free and open-source R environment (which is based on a proprietary programming language), as GEE is not designed for statistical research [38].
This study project falls within such a scientific framework.Its goal is to evaluate the new Landsat 9 data's capability for measuring albedo and, hence, investigate its utility in environmental monitoring.To do this, the algorithms that were proposed by da Silva et al. [30] and Liang [28] for estimating albedos from Landsat images were applied by writing two appropriate programs in the GEE and R environments.The approach was assessed using Landsat 8 and 9 images from two cities in Southern Italy: Palermo (Sicily), and Cagliari (Sardinia).Lastly, the impact of the Copernicus Urban Atlas 2018 land use/land cover (LU/LC) classes on data performance was investigated as well [10,[39][40][41].

Study Sites
The Mediterranean coastal cities of Palermo and Cagliari in Southern Italy (Fig. 1) were chosen as study areas because, for these sites, images with comparable acquisition data were identified for both the Landsat 8 and Landsat 9 missions.Additionally, these zones have comparable geographical, geomorphological, and climate characteristics, with hot and dry summers and mild winters [42,43].Palermo is located on the northern coast of Sicily Island on a densely urbanized wide plain surrounded by mountains to the south and west.Its elevation is about 150 m above sea level (a.s.l.) [44].Meanwhile, Cagliari is in the southern part of Sardinia Island and stretches across various hills (with an elevation of around 80-100 m a.s.l.)It is surrounded by wetlands (such as lagoons and ponds) to the west, north, and east [45].

Workflow
Figure 2 describes the approach that was adopted in this research.After our Landsat 8 and Landsat 9 TOA image selection, GEE was programmed with a JavaScript code to assist with cloud detection and masking.Following this, the albedos were determined for both the Palermo and Cagliari sites by applying da Silva et al.'s [30] and Liang's [28] algorithms.

Fig. 2. Operative workflow
As a consequence, the resulting albedo thematic maps were statistically examined and compared using R software [38], taking the entire Palermo and Cagliari sites as well as each LU/LC class derived from the Copernicus Urban Atlas 2018 into account [39].Tukey's filter [46,47] was then applied to detect and remove any outliers, which were explored and compared across both territories and each LU/LC class.

Input Data and Pre-processing Step
The input satellite data was chosen based on two key criteria: less-than-10% cloud coverage, and the same collecting day and hour.The first criterion assured that the clouds had the least influence on the processing output, while the second provided equivalent illumination conditions and, therefore, a similar performance comparison.This last criterion was satisfied since an underfed maneuver was completed between November 11 and 17, 2021 (the Landsat 9 satellite traveled beneath the Landsat 8 spacecraft [22]).Thus, the images that were captured on November 15 and November 13, 2021, for Palermo and Cagliari, respectively, matched both requirements and were selected.Their characteristics are described in Table 1.To determine the albedos, TOA images that had previously been orthorectified and calibrated were used.As a result, no orthorectification or radiometric calibration techniques were performed (other than cloud masking).In addition, the Palermo and Cagliari LU/LC Urban Atlas 2018 shapefiles [39] (Fig. 3) were imported in GEE.Urban Atlas (UA) is a local data-collection site that is provided by the Copernicus initiative, which is supported by the European Space Agency (ESA) and the European Environment Agency (EEA).The initiative features 28 LU/LC classes -17 of which are urban, and 11 of which are rural [48].The LU/LC class distributions (in percentages) are reported in Figure 4.

Land Surface Albedo Estimation
Da Silva et al.'s [30] and Liang's [28] algorithms were used to determine the broadband albedos from the Landsat 8 and various EO satellite data, respectively.The first approach was based on Equation (1), which was introduced by Zhong and Li [49] and Bastiaanssen et al. [50]: where α TOA , α ATM , and τ denote the TOA albedo, the path radiance albedo, and the atmospheric transmissivity, respectively.According to [51], α ATM was between 0.025 Source: [39] a) b) and 0.04; as recommended by [52], a value of 0.03 was used.The value of τ was, instead, determined by using Equation ( 2 where P o , W, and K t are the local atmospheric pressure [kPa], the precipitable water [mm], and the air turbidity coefficient, respectively.P o was measured at 10:00 on the days of the satellite images by the Palermo and Cagliari stations of the Italian National Tidegauge network (https://mareografico.it/).W was obtained from the "National Center for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) Reanalysis data, Water Vapor" data set [54] (available in the GEE catalog), which provides W values at six-hour temporal resolutions (00:00, 6:00, 12:00, and 18:00).The average of the two closest-in-time values of W (6:00 and 12:00) was computed to select a value of W that was representative of the Landsat acquisition time (around 10:00) for each satellite image (as described in [55]).A K t value of 1 is used for pure air, while a value of 0.5 is used for polluted air.Because the regional authorities detected air pollution levels at both sites that were lower than the legal limits in November 2021, a value of 1 was chosen in this case.Lastly, the Sun's zenith angle (Z) that was reported in the image metadata was employed.
The TOA albedo was, instead, computed by linearly combining the TOA reflectance (r i ) for each OLI-Landsat 8 band through the weighting coefficient (p i ) that is reported in Table 2  p r p r p r p r p r p r α = ⋅ + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ where p i is obtained by calculating the ratio between the solar constant of each spectral band (k i ) (computed using Equation ( 4)) and the sum of all of the solar constants that are used in the albedo computation (k TOT ) [10,30]: where L i and d are the radiance for each band and the correction of the eccentricity of the terrestrial orbit, respectively.On the other hand, Liang proposed Equation (5) to compute the albedos from the Landsat atmospherically corrected satellite images [28].Table 3 shows the weighting coefficient (ω i ), applied to the surface reflectance (ρ i ), that was adopted in this calculation:

Statistical Analysis
The resultant albedo maps were sampled and statistically investigated in the R environment in order to (i) evaluate the albedo distribution and variability within the study areas and LU/LC classes using base statistics metrics, (ii) investigate the Landsat 9 data potentialities in estimating the albedos using algorithms that were optimized for Landsat 8 and comparing their results to those that were obtained by processing the Landsat 8 images, and (iii) assess the outliers' impacts on the generated albedo maps.To meet the first objective, the mean (μ), standard deviations (SD), median (m), distribution trends, kurtosis (κ), and skewness (sk) were computed.For the second one, scatterplots that compared the Landsat 8 and Landsat 9 estimates were calculated, along with the correlation coefficient (ρ) and root mean square error (RMSE).Lastly, the third purpose was addressed by using Tukey's statistical approach; this allowed for an effective technique for reducing the inherent noise by locating and filtering out the outliers.According to Equation ( 6), an observation is labeled as an outlier when its value exceeds the lower (Q 1 ) or upper quartiles (Q 3 ) [46][47]: The above-mentioned statistical metrics were recalculated from the filtered albedo maps for the entire study areas and for each LU/LC class.These results were compared to those that were obtained while dealing with the unfiltered maps.

Land Surface Albedo Maps and Statistics Metrics
The albedo maps were produced to highlight the variability of the albedos within the case study territories.In Figures 5 and 6, the darker regions show lower values, whereas the yellow-orange areas have higher values.The lowest albedo values in Palermo could clearly be observed in the mountainous areas (mainly identified as the "Forest" and "Herbaceous vegetation" LU/LC categories), while in Cagliari, these were associated with the "Water" class.Additionally, the highest values could be found in the urbanized areas as well as in the classes such as "Mineral extraction and dump sites," "Open spaces," and "Herbaceous vegetation" for both sites.Table 4 presents the statistical metrics that were calculated for the entire territories of Palermo and Cagliari by using the two methods that were described in the previous sections.Based on the results of da Silva et al.'s technique, the μ value was around 0.2 in each map; moreover, the SD in Palermo was 0.074 and 0.071 for the Landsat 8 and Landsat 9 estimations, respectively.In Cagliari, the SD values for the Landsat 8 and Landsat 9 maps were 0.071 and 0.067, respectively.Furthermore, ρ was quite high in both areas (reaching 0.972 for Palermo and 0.942 for Cagliari), while the RMSE values were low (totaling 0.020 for Palermo and 0.024 for Cagliari).On the other hand, Liang's algorithm showed lower values than the other method; in fact, its μ value was approximately 0.1 in both study areas, and its corresponding SD was around 0.5.Nonetheless, the albedo values for Palermo were higher than those that were obtained in the Cagliari area according to da Silva et al.'s algorithm.
Histograms were created to determine how the data was distributed (Fig. 7).These distributions were visually congruent with the previously reported metrics values (Table 4).
Figure 8 depicts boxplots of the albedo values, which illustrate that many outliers were present in the upper section of the graph.Therefore, this signified that high albedo values (greater than 0.4 and 0.3 for da Silva et al.'s and Liang's algorithms, respectively) were detected for both data types and for both study areas.Instead, the scatterplots between the Landsat 8 and Landsat 9 results for both sites are illustrated in Figure 9.
Most of the spots in Palermo (Fig. 9a, c) were scattered along the diagram's bisector line; this indicated that many of the albedo values from Landsat 8 and Landsat 9 were linearly correlated.However, there were some inaccurate values -particularly in the lower-left corner of the image.Cagliari's pattern (Fig. 9b, d) was similar to Palermo's, but there were more mismatched spots across the plot.This is verified by Table 4, which indicates that ρ was lower in Cagliari than in Palermo and that RMSE was greater in Cagliari than in the corresponding Palermo area.The key statistical metrics for each LU/LC class in relation to the Palermo site are provided in Tables 5 and 6.The same data that was obtained for the Cagliari classes is shown in Tables 7 and 8. Furthermore, the scatterplots demonstrate the presence of multiple albedo values that were greater than 0.4 and 0.3, which corresponded to the outliers in the boxplots for da Silva et al.'s and Liang's approaches, respectively.To detect and remove these outliers as well as explain their impacts on the calculations, Tukey's filter and an outlier analysis were employed.According to Tables 5 and 6, Palermo had very high ρ values (>0.86) and low RMSE values (≤0.030) for all of the classes.On the other hand, Cagliari had four classes with ρ levels that were lower than 0.80 and/or RMSE levels that were greater than 0.03 (Tables 7, 8).

Impact of Tukey's Filter on Statistics Metrics
Because each study site had a high number of outliers (as shown in Figure 8), Tukey's filter [42] was applied to detect them.The statistical measures were then recalculated.As a consequence, Figures 10 and 11 and Table 9 show the distribution trend, scatterplots between the Landsat 8 and Landsat 9 images, and the statistical metrics for the whole region, respectively.Conversely, the metrics that were computed for each LU/LC classes are reported in Tables 10-13.Tukey's filter allowed us to raise the ρ values to greater than 0.952 and reduce the RMSE values to less than 0.019 for both areas and methods (Table 9).Furthermore, the skewness went negative in both sites when using Liang's technique and in Palermo when using da Silva et al.'s approach, while the kurtosis of the albedo distributions remained positive in both study regions (Figure 10 and Table 9).
The scatterplots of the filtered distributions indicate that the bulk of the miscorresponding albedo values were still preserved, but the higher values (greater than 0.4 and 0.3 for da Silva et al.'s and Liang's algorithms, respectively) were removed (Fig. 11).In line with the statistical metrics that were identified for the entire territories (Table 9), Tukey's filter increased the ρ value (greater than 0.88 and 0.86 for da Silva et al.'s and Liang's algorithms for both study sites, respectively) for the LU/LC classes (Tables 10-13), and it reduced the RMSE values (lower than 0.20 for Palermo, and 0.22 for both sites) for the LU/LC classes (Tables 10-13).

Outlier Analysis Output
Figures 12 and 13 represent the outliers in the Palermo and Cagliari regions that were obtained after using da Silva et al.'s and Liang's algorithms, respectively.
Most of the observed outliers in both areas overlapped (Figs.14,15).In Palermo, overlapping outliers represented more than 89% of the totals for both images and the applied method; in Cagliari, these accounted for more than 76% and 43% of the totals for da Silva et al.'s and Liang's methods, respectively.Palermo and Cagliari LU/LC class.The outliers in Palermo tended to be found in "Industrial, commercial, and public, military, and private units", "Herbaceous vegetation," and "Mineral extraction and dump sites," whereas the classes with most outliers in Cagliari were "Industrial, commercial, and public, military, and private units," and "Port areas" (Figs.[16][17][18][19].

Discussion
This research aims to test da Silva et al.'s and Liang's algorithms [30] for assessing the potentialities of the new Landsat 9 data for retrieving albedos in a Mediterranean urban context.For this purpose, the GEE platform and R open-source software were employed for the satellite-image processing and statistical analysis, respectively.
A first visual examination of the albedo maps (Figs. 5, 6) indicated that a comparable pattern of albedo values occurred across the territories of the investigated study sites (Palermo and Cagliari).Many pixels on the maps appeared to have albedo values that were between 0.1 and 0.3, although some discrepancies between the two regions were possible owing to differences in the geomorphologic and LU/LC features.Furthermore, the albedo-value distribution that was derived from Landsat 8 and Landsat 9 was consistent for both sites and for both of the applied algorithms.These assumptions were supported by the statistical analysis, which was valuable for retrieving the mean values and variability of the albedos over the two areas as well as for understanding the correspondences and discrepancies between the Landsat 8 and Landsat 9 estimations.
According to the statistical metrics that are reported in Table 4 for both case studies and algorithms, the Landsat 8 satellite data produced slightly higher μ, SD, and m values than the Landsat 9 data did (with a difference that was less than 10 −2 ).The only exception was Liang's algorithm for the Palermo study region; however, the albedo values that were obtained using Liang's technique were lower than those that were achieved using da Silva et al.'s approach according to the literature.Indeed, the mean albedo values that were determined using da Silva et al.'s and Liang's approaches for both Landsat 8 and Landsat 9 discoveries and the study areas were quite close to 0.20 and 0.10, respectively (Table 4).This was most likely due to the two cities' comparable climates and latitudes.Furthermore, this conclusion corresponded to those that were presented in previous research publications when combined with the SD values: Taha [56] indicated that the albedo values in urban regions ranged between 0.10 and 0.20 for several European cities, while Brivio et al. [57] reported that the albedo values in urbanized areas ranged between 0.12 and 0.21.Secondly, the ρ values were very high, whereas the RMSE values were low for both study sites and the applied algorithms (Table 4).This is supported by the scatterplots in Figure 9, which indicate an excellent correlation of the albedos for both regions (despite the fact that Cagliari had a higher number of mis-corresponding values).As a result, the Palermo region had a greater correlation and a smaller RMSE between the Landsat 8 and Landsat 9 estimations.
A combined visual inspection of the albedo maps and satellite images revealed that the mis-corresponding values for both areas were primarily due to the presence of small clouds (probably vapors) that did not coincide in the Landsat 8 and Landsat 9 images; alternatively, these were due to the radiometric differences in some pixels between the Landsat 8 and Landsat 9 images.According to Gross et al. [22] (who performed the cross-calibration of the Landsat 8 and Landsat 9 instruments using images that were collected during the underfly event in November 2021), the three main sources of uncertainty were i) geometric (due to the pointing differences between the two sensors), ii) spectral (due to the spectral alterations between the two sensors), and iii) angular (caused by the differences in the viewing/illumination angles).Furthermore, Cagliari had more small clouds (vapors) than Palermo, which might explain why the Cagliari scatterplot had a higher number of mis -corresponding albedo values (Fig. 9).As a result of the statistical comparison, da Silva et al.'s method [30] could be applied to Landsat 9 images with a good approximation -even if some radiometric differences were detected between them.This conclusion can be extended to Liang's method as well.Despite the fact that the latter method understated the albedo values, a strong correlation was found for both research locations between the Landsat 8 and 9 pictures.
Tables 5 and 6 highlight that the LU/LC classes for the Palermo site did not affect the main statistical metrics.In contrast, some issues were detected in Cagliari, where the RMSE value was greater than 0.030 for several LU/LC classes (such as "Industrial, commercial, and public, military, and private units," "Arable land," "Discontinuous urban fabric," and "Wetlands," for instance) (Tables 7, 8).Similarly, the ρ values were below 0.065 with the employed approach for three categories ("Arable land," "Discontinuous low density urban fabric," and "Wetlands").The "Wetlands" category (which was not present in the Palermo area) was the most critical one.Moreover, the Landsat 8 data revealed relatively higher μ values than the Landsat 9 data for all of the LU/LC classes in Palermo and Cagliari; this was in accordance with the statistical metrics that were associated to the entire territories.In terms of SD (and solely for the Palermo area), the Landsat 8 satellite images produced somewhat higher values than the Landsat 9 data for most of the LU/LC categories.In Cagliari, however, the Landsat 9 SD values were found to be greater than those that were derived from the Landsat 8 data for numerous classes.
Previous studies have provided typical albedo values for the various LU/LC categories.According to Allen et al. [51], for instance, agricultural landscapes had albedo values that ranged from 0.14 to 0.22, and water bodies had albedo values that ranged from 0.025 to 0.348.According to Brivio et al. [57], roads had albedo values that ranged between 0.10 and 0.28, bare soils -between 0.05 and 0.31, and vegetation -between 0.14 and 0.45 (depending on its phenological stage).When combined with the SD values, the main albedo values that were obtained in this study for the various classes in both the Palermo and Cagliari sites (Tables 5-8) were comparable to those that were reported in the aforementioned contributions.Figure 8's boxplots reveal that there were several outliers in both of the Palermo and Cagliari study areas; these outliers were identified and removed using Tukey's filter.As a result, the original albedo datasets were filtered, and all of the statistics were recalculated.Even after Tukey's filter, the Landsat 8 and Landsat 9 outputs had roughly identical values of μ, SD, and m (with somewhat higher values being produced from the Landsat 8 data) (Table 9).Nonetheless, ρ was significantly improved and RMSE slightly reduced for both of the used methods for Cagliari.For Palermo, the RMSE and ρ values that were derived from da Silva et al.'s technique were a bit lower and slightly greater, respectively, than those that were obtained from the unfiltered pictures.In fact, Tukey's filter excluded only those outliers with albedo values that were greater than 0.4; therefore, Palermo (which had mis-corresponding values that ranged from 0 to 0.25 -Fig.9) partially benefited from its application.Furthermore, the kurtosis decreased in both study areas after implementing the filter (albeit still positive); however, the skewness became negative in Palermo for da Silva et al.'s algorithm and both Palermo and Cagliari for Liang's algorithm (Fig. 10).The filtered distribution scatterplots (Fig. 11) demonstrate that Cagliari has greater mis-corresponding values between the Landsat 8 and Landsat 9 estimations than Palermo did; however, the difference was less in terms of ρ.As shown in Tables 10-13, Tukey's filter results indicated a minor drop in RMSE and an increase for the majority of the analyzed LU/LC categories in both Palermo and Cagliari.
As detailed in Chapter 3, the positions of observed outliers were determined and mapped.Because most of them overlapped in both study zones (Figs.14, 15), a connection between the Landsat 8 and Landsat 9 estimations was not reached for a tiny percentage of them.With the exception of a 2-3% variation, the performances of the two satellite sensors may be considered to be similar.Nevertheless, Liang's method detected more non-overlapping outliers than the other approach did (Fig. 15).However, similar results were extracted from the Landsat 8 and 9 images.
An outlier study based on the investigated LU/LC categories further confirmed these outcomes (Figs.16,17).In Palermo, outliers from both the Landsat 8 and Landsat 9 data were largely identified in the "Industrial, commercial, and public, military, and private units," "Herbaceous vegetation," and "Mineral extraction and dump sites" categories.The first two of the mentioned classes were also among the most widespread in Palermo (Fig. 3a).On the other hand, Figures 18 and 19 demonstrate that, in Cagliari, the outliers were mostly located in the "Industrial, commercial, and public, military, and private units" and "Port areas" classes.Even in this scenario, the results remained the same when both the Landsat 8 and Landsat 9 images were used.As a result, the fraction of the outliers that were obtained from the Landsat 8 and Landsat 9 data (as well as the findings for both study regions) were almost comparable.

Conclusion
The purpose of this research was to evaluate the accuracy of the Landsat 9 data for albedo estimation in urban contexts.Da Silva et al.'s approach [30] (created for the Landsat 8 data) and Liang's algorithm [28] (introduced for handling various forms of EO satellite data) were evaluated in order to achieve this study objective.
Palermo and Cagliari (cities in southern Italy) were selected as test sites because it was possible to acquire Landsat 8 and Landsat 9 images on the same day and within a suitable time frame.Such open medium-resolution satellite data is preferable for environmental monitoring and climate change studies, as it allows us to estimate many bio-geophysical variables as well as their spatial and temporal variabilities [10,16,[57][58][59][60][61][62][63][64][65].Measuring and monitoring albedo levels is crucial for addressing issues such as UHI and drought, which are caused by changes in the energy fluxes between land surfaces and the atmosphere [1,3,5,[10][11]18].
A comparison of the statistical measures that were produced from the Landsat 8 and Landsat 9 albedo maps showed an incredibly good correlation and very low RMSE for both study sites and the examined technique.With a few exceptions (most likely due to radiometric discrepancies or the existence of some tiny clouds [possibly vapors] in the Palermo and Cagliari areas), the results were largely encouraging.Such discrepancies resulted in a few mis-matched albedo values between the Landsat 8 and Landsat 9 computations.As a consequence, da Silva et al.'s [30] and Liang's approaches appeared to be significantly relevant to the Landsat 9 data; however, it should be highlighted that substantial radiometric variations existed for some pixels between the images that were taken by the two platforms.This work does not intend to lay out the best technique for extracting albedo values from Landsat 8 and 9 data but rather to establish the potential of using the same algorithms on both images.The comparability of their outputs implies that such tried-and-true procedures may be extended to Landsat 9 and that the data that is generated by this satellite platform can be utilized in the Landsat 8 continuity.
Lastly, the GEE cloud-based platform has been demonstrated to offer advantages over alternative desktop software when working with large amounts of geospatial data, as it decreases the acquisition and operating times.Thus, it appears to be the ideal tool for addressing major environmental challenges through the processing and analysis of big remotely sensed data.

Fig. 1 .
Fig. 1.Palermo (a) and Cagliari (b) study areas, along with corresponding geographical framework Source: Google Earth for geographical framework

Fig. 5 .
Fig. 5. Albedo maps retrieved from Palermo Landsat 8 (a) and Landsat 9 (b) images as well as from Cagliari Landsat 8 (c) and Landsat 9 (d) images by applying da Silva et al.'s algorithm

Fig. 7 .Fig. 8 .
Fig. 7. Distribution trends of albedo values generated from Palermo Landsat 8 (a) and Landsat 9 (b) data and Cagliari Landsat 8 (c) and Landsat 9 (d) data using da Silva et al.'s algorithm as well as Palermo Landsat 8 (e) and Landsat 9 (f) data and Cagliari Landsat 8 (g) and Landsat 9 (h) data using Liang's algorithm

Fig. 11 .
Fig. 11.Scatterplots between albedos estimated from Landsat 8 and Landsat 9 images for Palermo (a) and Cagliari (b) using da Silva et al.'s algorithm and estimated from Landsat 8 and Landsat 9 images for Palermo (a) and Cagliari (b) using da Silva's approach:all findings were obtained after applying Tukey's filter (L8 -Landsat 8; L9 -Landsat 9)

Fig. 12 .
Fig. 12. Outlier positions in Palermo (a) and Cagliari (b) obtained using da Silva et al.'s approach Source: Google Earth

Table 1 .
Selected image features

Table 2 .
Weighting coefficient values of da Silva et al.'s algorithm

Table 3 .
Weighting coefficient values of Liang's algorithm

Table 12 .
Main statistical metrics obtained using da Silva et al.'s approach for each LU/LC class of Cagliari after applying Tukey's filter (L8 -Landsat 8; L9 -Landsat 9)