ISLAND: Interpolating Land Surface Temperature using land cover (2024)

3.1 Landsat 8 Data

We collect Landsat 8 Collection 2 Level-2 LST products from Google EarthEngine(Gorelick etal.,, 2017).Landsat 8 provides LST products at $30\text{\,}\mathrm{m}$ resolution and a revisit cycle of 16 days.Landsat 8 collects top-of-atmosphere (TOA) spectral radiance values via its band 10 Thermal Infrared Sensor (TIRS).LST is derived from TOA spectral radiance using the radiative transfer equation(RTE)-based single-channel algorithm(Malakar etal.,, 2018).ASTER GED(Hulley etal.,, 2015) was used to correct the effect of surfaceemissivity.Atmospheric effects were compensated by the Goddard Earth Observing System, Version 5 (GEOS-5) reanalysis data using the radiative transfer model MODTRAN 5.2.The RMSE of Landsat LST products is approximately$2.2\text{\,}\mathrm{K}$(Malakar etal.,, 2018).Please refer toMalakar etal., (2018) for a more detailed descriptionof the Landsat LST retrieval algorithm.We denote the input Landsat LST image as $\tilde{T}\in\mathbb{R}^{H\times W}$, where $H$ and $W$ denote the height and width of the LST image, respectively.

The Landsat LST product includes cloud mask information(CFMask)(Zhu and Woodco*ck,, 2012) on a per-pixel basis.CFMask provides pixel quality attributes for each Landsat 8 LST image, indicating the presence of cloud contamination.We build a pixel-wise binary occlusion mask $O\in\{0,1\}^{H\times W}$ for each Landsat LST image, where the pixel value is set to True if CFMask indicates that there is cloud, cloud shadow, or cirrus.From the occlusion mask, we calculate the occlusion factor $\theta\in[0,1]$ that measures the fraction of pixels occluded.Formally, we have $\theta=\frac{\sum_{\boldsymbol{p}}O_{\boldsymbol{p}}}{H\cdot W}$, where $O_{p}$ denotes a pixel of $O$.$\theta$ is an important metric that measures the severity of cloud contamination.A higher $\theta$ indicates more severe occlusion.

3.2 NLCD Land Cover data

For land cover data, we use the National Land Cover Database(NLCD)(Yang etal.,, 2018) from the U.S. Geological Survey (USGS).NLCD provides spatially referenced descriptive data on the characteristics of the land surface using a set of thematic classes (e.g.,urban, forest, and agriculture).NLCD is available for the CONUS region at $30\text{\,}\mathrm{m}$ resolution with a release cycle of once per 3years.We denote the NLCD land cover image as $L\in\mathbb{Z}^{H\times W}$.

3.3 In situ data

The in situ LST is obtained by inverting the upwelling component of Eq.(1)

4.1 Spatial Channel

In our approach, the prediction of occluded pixel values is informed by leveraging the information embedded in their surroundings.We term this approach the spatial channel.The intuition behind the spatial channel is that nearby objects that are in the same land cover class are likely to exhibit similar thermal properties.This approach filters the data based on land cover class, since the distribution of LST is dependent on land cover, as depicted in Fig.2.The computation of the spatial channel bears a close resemblance to bilateralfiltering(Paris etal.,, 2009).While traditional bilateral filtering considers variations in pixel intensities of the input image with the aim of preserving sharp edges, the spatial channel uses a different approach, where we filter the input image using the pixel attributes informed by auxiliary images (land cover class $L$ and occlusion $O$).

Consider our goal to be predicting $\tilde{T}$ at an occluded pixel location $\boldsymbol{p}$.Under a low occlusion factor $\theta$, we estimate $\tilde{T}_{\boldsymbol{p}}$ using a weighted average estimate from other pixels.Conceptually, the weights are determined based on three factors: proximity, land cover class label, and cloud occlusions.

Algorithm 1 shows the implementation of the spatial channel.Let $\tilde{T}$ be the occluded LST image as input, $f$ be the size of the local window, and $\theta^{*}$ be the maximum occlusion factor threshold for local filtering.We obtain interpolated image $\hat{T}_{sp}$ as follows:

In our implementation, we set the values of the optimization parameters as $f=75$ and $\theta^{*}=0.5$. These values were selected based on extensive testing on a subset of images in our dataset to achieve the best performance.

4.2 Temporal Channel

In this section, we show how to generate an interpolated prediction using temporal information.We call this method the temporal channel, which is complementary to the spatial channel defined in Sec.4.1.As seen in Fig.3, the intuition behind the temporal channel is that objects in the same land cover class tend to exhibit similar thermal dynamics over time.The temporal channel involves four steps: (a) select a set of frames as a reference, (b) preprocess each reference frame via the spatial channel, (c) apply linear adjustments to each reference image, and (d) interpolate occluded regions based on the set of adjusted reference frames.

Reference frame selection:We select reference frames based on two conditions: (i) seasonality and (ii) cloud occlusion.The goal is to identify suitable reference images that can be used to accurately reconstruct occluded LST.

For seasonality, we take into account the temporal offset between the occluded target image $\tilde{T}$ and other available images $\tilde{T}^{\prime}$.Due to temporal continuity, the previous and the next LST sample over the same region are good reference points.Constrained by the revisit cycle of Landsat 8, the previous and next samples are taken 16 days before and after the target date.In addition to the immediate temporal neighbors, we also leverage the predictable seasonal variations in LST signals. This allows us to include the previous and the next sample acquired in other years as potential references.In our actual implementation, we increase the selection limit to samples collected 2 cycles prior to or after the target date in each year, defining the temporal bracket duration $\delta$ as 2.See vertical gray stripes in Fig.3(b) for a visual illustration.

The second condition of reference frame selection is the occlusion factor $\theta$.We have observed that minimally occluded reference frames lead to smaller interpolation errors.Therefore we only consider selecting reference frames whose $\theta$ is below a certain maximum tolerable threshold, $\theta_{\text{max}}$ (see red dotted line in Fig.3(b)).We call the set of images satisfying these two conditions candidate reference frames $\mathcal{T}^{\ast}_{ref}$.Within $\mathcal{T}^{\ast}_{ref}$, we prioritize selecting images that are captured closer to the target frame in time, as measured by the temporal difference $\Delta(|\tilde{T}^{\prime}-\tilde{T}|)$.To achieve this, we choose a subset of $n$ frames from $\mathcal{T}^{\ast}_{ref}$ with the lowest $\Delta(|\tilde{T}^{\prime}-\tilde{T}|)$.We call this selected subset of reference frames $\mathcal{T}_{ref}$.

Spatial channel pre-processing:Given a selection of reference frames $\mathcal{T}_{ref}$, we apply the computation of spatial channel from Fig.4.1 to produce a set of spatially complete reference frames.As a result, each pixel contains either observed or interpolated temperature data.Note that the imposed constraint on $\theta_{max}$ means that the spatial channel only needs to interpolate a minimal amount of occlusion.

Linear adjustments:After pre-processing $\mathcal{T}_{ref}$, we apply a linear adjustment to all pixels in each class. Fig.3(a) shows that the changes in the mean LST (denoted as $\mu_{c}$) can differ drastically between classes.Although the selection of reference frames helps mitigate these discrepancies, some adjustments are still necessary.Specifically, for each reference frame, we add the difference in $\mu_{c}$ (i.e.,$\Delta\mu_{c}$) between two dates to all pixels belonging to the corresponding land cover class.$\Delta\mu_{c}$ visually translates to the length of each bar in Fig.3(a).We repeat this process for all images in $\mathcal{T}_{ref}$.

Reference frame-based interpolation:After the previous two steps, we obtain a set of linearly-adjusted reference frames.The interpolated LST image, $\hat{T}^{temp}$, is the average of each linearly adjusted reference frame.This interpolation step combines the information from multiple reference frames to produce a spatially complete and temporally consistent estimate of the occluded LST image.

Algorithm 2 shows the implementation of the procedures above. Let $n$ be the number of reference frames to be selected, $\delta$ be the temporal bracket duration, and $\theta_{max}$ be the maximum tolerable occlusion factor. We compute the interpolated image $\hat{T}_{temp}$ as follows:

We set $n=3$, $\delta=2$, and $\theta_{\text{max}}=0.1$. SP($\cdot$) denotes the function for spatial channel defined in Algorithm 1, and $f$ and $\theta_{\text{min}}$ at line 4 follow their default values from Sec.4.1.The subscript $\cdot_{(n)}$ used in line 2 signifies the first $n$ elements of the sequence, arranged according to the minimization criterion.$P_{c}$ is the set of all pixel locations where the corresponding land cover class is $c$, and $|P_{c}|$ is the cardinality of set $P_{c}$.

4.3 Estimate LST via Weighted Average

5.1 Public Data Products

We deploy ISLAND across 20 regions in the United States. These regions areselected for having the highest populations, as reported in the 2020 U.S.Census(United States Census Bureau,, 2023).³³3For each region, we manually define apolygon roughly covering the metro region for each city. Due to the 32 MB limitper image download from Earth Engine, the polygon may not encompass the entiremetro region for some cities, such as New York.We collect Landsat data from 2019–2023 and NLCD 2021 release using the data compilation process described in Sec.3.We then use ISLAND to predict cloud-contaminated LST for each region.ISLAND produces interpolated LST at a spatial resolution of $30\text{\,}\mathrm{m}$ every 16 days⁴⁴4Subject to data availability and requires $\theta<0.99$. for each observation region.Our public dataset is available on the NHERI DesignSafeCyberinfrastructure⁵⁵5https://doi.org/10.17603/ds2-3rf5-sd58(Rathje etal.,, 2017).

Fig.4 shows a set of interpolated LST maps produced by ISLAND. We choose a diverse set of examples consisting of different land cover and cloud occlusion conditions. ISLAND reconstructs spatially complete LST maps under a variety of conditions, from lightly occluded by thin cirrus clouds (Houston) to heavily occluded by optically thick clouds (New York, 88% occluded). Our model performs well across regions with different land cover characteristics, from dense urban settings (New York and Los Angeles) to diverse wetlands (Jacksonville).

5.2 Simulation Evaluation

To evaluate the performance of our LST reconstruction method, we simulate cloud contamination by artificially occluding Landsat LST images. The added occlusions occupy a set of $s\times s$ rectangular regions. We set LST pixel values to zero and occlusion mask to True in these regions. To evaluate, we choose four urban regions with different climatological and surface land cover conditions, as seen in Table1. For each region, we apply up to $n$ occlusion regions of size $s\times s$ pixels for all available LST images from 2019–2023. We use NLCD 2021 release as input for our simulation evaluation. For LST images with actual cloud occlusion (i.e.,real occlusion), we place artificial occlusions alongside real occlusions and ensure no overlap. The evaluation metrics are computed only in artificially occluded areas by comparing them to the original LST images. We report the mean absolute error (MAE) and root mean squared error (RMSE) in Table1. M1 refers to our model ISLAND, while M2 - M4 refers to other models discussed in 5.3.

Table1 suggests that the RMSE typically ranges from 1.65–$2.62\text{\,}\mathrm{K}$ for a variety of urban regions and occlusion scenarios. The first two rows of Table1 indicate that more occlusions generally lead to larger reconstruction errors. ISLAND performs well in dense urban regions (such as New York City) and diverse wetlands (such as Jacksonville). Our simulation evaluation indicates that ISLAND is robust in performance and is able to generalize reasonably well to different regions in the U.S. and under different occlusion characteristics.

5.3 Ablation Study on Simulation Data

To further demonstrate the effectiveness of our model, we perform an ablation study⁶⁶6Not to be confused with the glaciological definition of ablation, which refers to the process of removing snow, ice, or water from a glacier or a snow field. Here, we use the artificial intelligence definition of ablation study, where certain components of a model are removed in order to gain a better understanding of the model’s behavior. , where we remove key components of our model and observe the impact of each of these components on the overall performance. Table1 shows a list of models. M1 refers to our full model, ISLAND. In M2, we exclude the temporal channel and only keep the spatial channel. In M3, we discard the spatial channel and only keep the temporal channel. Note that M1 is a weighted average of M2 and M3, following Eq.(8). M4 is the same as ISLAND, except we remove NLCD land cover labels as input. M5 adopts a simplified approach where missing values are replaced with the average pixel value within an image without considering their spatial distribution or land cover class.

Table1 shows that the use of NLCD land cover labels (M1) leads to better reconstruction performance over the model without NLCD as input (M4). Such difference is more pronounced in dense urban regions, such as New York, and heterogeneous regions, such as Jacksonville.

When confronted with heavier occlusions (large $\theta$), M3 consistently emerges as the top-performing model.This finding highlights the effectiveness of relying solely on the temporal channel in such challenging scenarios.Conversely, for occlusions of low to moderate severity, the spatial channel exhibits stronger performance.By appropriately favoring the temporal prediction in the presence of heavy occlusions, and the spatial prediction for lighter occlusions, our weighting scheme (Eq.(8)) allows for adaptability to varying occlusion levels, enhancing the model’s robustness and accuracy in capturing LST under diverse conditions.

5.4 In situ Evaluation

In this section, we evaluate ISLAND reconstruction results against in situ LST measurements collected in four SURFRAD stations. Our in situ evaluation uses Landsat and SURFRAD data collected in 2013–2020 and NLCD 2016 release as inputs to ISLAND. In 5.2, simulation evaluation examines ISLAND’s ability to reconstruct a theoretical clear-sky LST under artificial occlusion. However, studies have shown that clouds have a cooling effect on the shaded region(Weng and Fu,, 2014), and simulation data is unable to evaluate if our reconstruction algorithm accounts for this effect. As such, we use in situ measurements to evaluate reconstruction performance under cloud-sky conditions.

There are two primary disadvantages to evaluating ISLAND against SURFRAD in situ data.First, remote retrieval of LST, even in clear-sky conditions, has uncertainties.Sources of error include uncertainty in emissivity estimation, atmospheric compensation, etc.Further, a Landsat 8 pixel (at $30\text{\,}\mathrm{m}$ resolution) is larger than thefield-of-view (FoV) of the SURFRAD instrument(Malakar etal.,, 2018).Differences in FoV compounded with spatial heterogeneity in temperature atSURFRAD sites(Malakar etal.,, 2018) cause Landsat LST to furtherdeviate from in situ LST. Therefore, we report Landsat LST versus in situ LST under clear-sky conditions (without involving ISLAND) in Fig.5. RMSE ranges from $2.95-$$2.99\text{\,}\mathrm{K}$ across four sites.⁷⁷7We choose not to include two other SURFRAD sites, BND and FPK, because the Landsat clear-sky RMSE at these sites is too high in our calculation, at $3.55\text{\,}\mathrm{K}$ and $4.78\text{\,}\mathrm{K}$ respectively. The error here represents the underlying uncertainties of comparing Landsat LST with in situ LST, which is external to our reconstruction algorithm.

The second disadvantage of SURFRAD evaluation is that all SURFRAD sites are located in rural, hom*ogeneous areas. Recall that the primary advantage of ISLAND is the use of NLCD data, which is applicable to urban areas but not to SURFRAD sites. Unfortunately, there are no publicly available in situ LST validation sites located in urban regions. Despite these limitations, we report in situ validation results to build an understanding of how ISLAND performs under real cloud occlusion.

Fig.6 shows ISLAND reconstructed LST versus in situ LST under cloudy-sky conditions. We define a given data point as cloudy-sky if and only if the corresponding Landsat pixel is flagged as cloud, cloud shadow, or cirrus, according to CFMask. Reconstruction RMSE ranges from $3.90-$$4.75\text{\,}\mathrm{K}$, across four sites. Compared to clear-sky RMSE in Fig.5, ISLAND introduces an additional $0.95-$$1.76\text{\,}\mathrm{K}$ RMSE error across four sites, with the average additional RMSE being $1.33\text{\,}\mathrm{K}$.

Fig.6 also shows that there is no systematic overestimation (i.e.,positive bias) across all four SURFRAD sites, suggesting that ISLAND effectively accounts for the local cooling effects caused by clouds. We believe that the modeling of the local cooling effect is primarily driven by the spatial channel defined in 4.1. As clouds move, they also cool the surrounding neighborhood. Under the assumption that surface objects have some degree of thermal inertia, we believe that the spatial channel utilizes surrounding cooler pixels to account for the local cooling effect caused by clouds, thereby accurately predicting LST under cloudy-sky conditions. While simulation evaluation (Table1) shows that using the temporal channel alone (M3) leads to better results than M1, in situ evaluation, however, suggests that using both spatial and temporal channels (M1, Fig.6) leads to better results. When using temporal channel only (M3), the RMSE at DRA, SXF, PSU, and GWN are $3.96\text{\,}\mathrm{K}$, $5.07\text{\,}\mathrm{K}$, $4.88\text{\,}\mathrm{K}$, and $3.86\text{\,}\mathrm{K}$, respectively.

Finally, Fig.6 includes reconstruction error for all cloudy-sky conditions except for dates with more than 99% of pixels occluded ($\theta>0.99$), demonstrating ISLAND’s robustness under a wide variety of cloud occlusion scenarios, including severely occluded LST images.

5.5 Applications

In this section, we show a set of applications demonstrating the impact of ISLAND on a variety of LST applications in urban environments.

5.5.1 Deriving surface temperature trends

The robustness demonstrated in Sec.5.2 – 5.4 and the compelling results displayed in Fig.4 underscore the model’s ability to generate accurate interpolated LST values across a wide range of conditions, with the only constraint being that the occlusion factor $\theta<.99$. As highlighted in Sec.1, previous studies investigating changes in land surface temperature through remote sensing have been constrained by limited observational conditions, restricting their analyses to a fraction of dates characterized by minimal cloud occlusion (Sobrino etal.,, 2013; Baiocchi etal.,, 2017; Huang and Wang,, 2019; Gomez-Martinez etal.,, 2021). However, with the introduction of ISLAND, these limitations become obsolete, granting access to a significantly expanded set of operational LST data, particularly in urbanized regions.

Beyond the production of interpolated image outputs, ISLAND enables the examination of temporal variations in LST on daily-to-seasonal timescales. By reconstructing skillful LST maps for the majority ($\theta<.99$) of observation dates, ISLAND enables the comparison of thermal behaviors for a given region across time, at a relatively dense sampling rate of every 16 days, and encompassing diverse land cover types. To illustrate this capability, Fig.7 showcases the evolution of surface temperature in Houston. Each colored line represents the average LST for a particular land cover class for all grid cells over Houston. Pronounced temperature seasonality is evident in the time series from 2019–2023, with clear changes in seasonal structure from year to year. The ability to partition surface temperatures retrieved from different land cover types reveals differences of up to $15^{\circ}$C between forested, water-covered surfaces and urban developed surfaces (e.g.,open-water vs. developed high intensity). With a 16-day temporal resolution, one can evaluate temperature distributions and variances over different land cover classes in different seasons. ISLAND facilitates an in-depth investigation of the temporal dynamics of surface temperature within any region located in the CONUS, providing valuable insights for climatological and ecological analyses.

5.5.3 Social Vulnerability & Urban Heat Stress

As illustrated in the last two applications, ISLAND facilitates the development of comprehensive datasets of LST and UHIE. The developed datasets will enable better characterization of heat exposure and its impacts on social, infrastructure, and environmental systems. A representative example application would be to investigate inequities in urban heat exposure.Given the health, well-being, and quality of life implications of urban heat,and initiatives like Justice 40(The White House,, 2022), whichcall for federal climate investments to be directed to environmental justicecommunities, understanding the equities in urban heat exposure can centrallyguide prospective investments.For example, quantifying inequities in exposure to urban heat will help design adaptation measures such as increasing vegetation cover or guiding urban planning, among others.

The distribution of UHIE for residential areas and the social vulnerability ofthe exposed population for a few cities are shown in Fig.9. Here,social vulnerability is measured using the Centers for Disease Control andPrevention Social Vulnerability Index (CDC SVI)(Centers for Disease Control etal.,, 2020). TheCDC SVI measures social vulnerability on a scale of 0 (least vulnerable) to 1(most vulnerable), taking into account socioeconomic status (e.g.,housing costburden), household characteristics (e.g.,civilian with a disability), racialand ethnic minority status (e.g.,Hispanic, Alaska Native), housing type andtransportation factors (e.g.,no vehicle). The latest available residentialland use data (2016) fromMcShane etal., (2022) are used to identifyresidential regions. Additionally, UHIE is calculated as the number of daysthat a pixel (resolution of $30\text{\,}\mathrm{m}$ $\times$ $30\text{\,}\mathrm{m}$)exceeds a land surface temperature threshold of $35^{\circ}$C(308.5K). Only pixels with at least one day of temperature over the threshold areconsidered for the analysis.

From Fig.9, many cities show a systemic inequity in heat exposure. For example, in Los Angeles, socially vulnerable communities are exposed to high urban temperatures compared to less socially vulnerable communities (Pearson’s correlation, $r=0.58$). A similar trend can be seen in cities such as San Antonio ($r=0.45$) and San Francisco ($r=0.46$), San Jose ($r=0.44$), and New York ($r=0.34$). In contrast, cities such as Houston ($r=-0.01$), Dallas ($r=-0.10$), Jacksonville ($r=0.17$), and Fort Worth ($r=0.08$) show no or negligible inequity in urban heat exposure. Factors ranging from vegetative cover to colocation with industrial or commercial locations might influence urban heat exposure. By providing reliable and more complete datasets to quantify UHIE, this study will allow for a better understanding of the factors influencing heat exposure and, as a result, will aid in developing strategies to mitigate urban heat stress.

While this section only provides three case studies, the proposed method facilitates many applications requiring high-resolution site-specific data. Some example applications include designing building envelopes and Heating, Ventilation, and Air Conditioning (HVAC) systems,investigating the influence of thermal stresses on infrastructure aging and deterioration, power system reliability, and energy demand shifts,or even high-resolution weather and natural hazard modeling that accounts for fine-scale surface temperature effects.

6.1 Assumptions and Limitations

Despite the demonstrated validity and effectiveness of ISLAND, the model does contain important assumptions and limitations.

Limitations of NLCD land cover labels:As stated in Fig.1, our model uses the NLCD land cover data.Specifically, we use the NLCD 2021 release(Dewitz,, 2021) for ourpublic data products, which most accurately reflect the state of land coverlabels in the year 2021.Our Landsat LST inputs, in contrast, span from 2019–2023.Within this observational period, NLCD 2021 alone does not reflect changes in land cover due to urban expansion, meandering, or coastal erosion.NLCD releases every three years, and there are versions available for 2019, 2016, etc.Although a fraction of observation dates could have potentially benefited from using the 2019 release instead of 2021, we chose not to implement computation using multiple NLCD versions for simplicity.Despite using only the NLCD release for 2021, we still observe a significant advantage in terms of reconstructed LST when employing the NLCD dataset as input, as seen inTable1.

Reliance of accurate cloud masks:Another assumption of our model is that all clouds are correctly labeled. Our algorithm relies on accurate cloud labeling, and errors in cloud labeling would most likely propagate to LST estimations.Recall from Fig.1 that we utilizeCFMask(Zhu and Woodco*ck,, 2012) to determine if a given pixel is affected bycloud, cloud shadow, or cirrus.CFMask algorithm can produce inaccurate cloud labels, leading to erroneous LST reconstructed values.Cirrus is a category of clouds known to be challenging todetect(Qiu etal.,, 2020).Row two of Fig.4 shows that unidentified cirrus in the lower left corner leads to visible discontinuities in the reconstructed LST image.Unidentified cirrus pixels seem to affect the spatial channel more than the temporal channel.Moreover, it is difficult to reliably detect clouds in snow-covered terrain dueto the spectral similarity between cloud andsnow(Stillinger etal.,, 2019).Consequentially, in a simulated evaluation, we observe significantly higher RMSE in Denver, a region with prolonged snow coverage.Finally, in some cases, we also observe that optically bright and thermally cold buildings are falsely labeled as clouds, though this mislabeling is rare.The increasing prevalence of white roofing applied to increase urban albedo anddecrease the UHIE may worsen the future uses of ISLAND(Fayad etal.,, 2021).It is important to address these challenges and improve cloud labeling techniques to ensure the accuracy of our model’s predictions.

Errors external to interpolation:As mentioned in Sec.5.4, remote LST retrieval itself has underlying uncertainties, even in clear-sky conditions.Other sources of error include, but are not limited to, sensor calibration,atmospheric profiles, and emissivityestimation(Li etal.,, 2013).Profiling and mitigating these types of errors are areas of activeresearch(Li etal.,, 2013) in the field of remote sensing; suchimprovements are beyond the scope of this paper.

6.2 Transferability and Scalability

We designed ISLAND to be easily accessible to the broader research community.All required inputs of our model listed in Fig.1 areacquired from publicly available sources and are extracted from Google EarthEngine(Gorelick etal.,, 2017) and its Python API, geemap(Wu,, 2020).Since our model is essentially based on a set of filters, we do not require massive computing power or GPU acceleration.For context, it takes roughly 2min to process one imageon a 12-core CPU (AMD Ryzen 5900) with 32GB of DRAM.

Currently, our model is only available to the CONUS region, constrained by the NLCD dataset.In Fig.6.3, we discuss potential avenues for expanding the operational region.

6.3 Potential Improvements and Future Opportunities for ISLAND

Extending area of study:In this paper, visual results (Fig.4), simulation evaluations (Table1), and demonstrated applications are all focused on urban regions.While in situ evaluation (Fig.5 and 6) are conducted on SURFRAD sites, which offers an indicator of model performance on regions with relatively high land cover hom*ogeneity (see Fig.A for visual illustrations), extensive future testing is required to better quantify performance outside of urban settings.

In Sec.6.1, we showed that the use of the NLCD dataset currently restricts our analysis to the CONUS region.In theory, we can potentially expand to other regions where there are appropriate land cover labels.For example, the Copernicus CORINE Land Coverdataset(Buchhorn etal.,, 2020) wouldfacilitate the extension of the model over Europe; the China land cover dataset(CLCD)(Yang and Huang,, 2021) is also available, enabling research over theAsian continent. It is important to note that the performance of ISLAND isimpacted by spatiotemporal resolution, diversity, and accuracy of land coverlabels.Therefore, careful consideration should be given to the suitability and quality of the land cover datasets when expanding the study area beyond CONUS.Additional testing is required to quantify ISLAND performance across different datasets and cities outside of CONUS with different land cover characteristics.

Improving temporal resolution:Higher temporal resolution for LST products is instrumental for downstream operational studies.Our model uses Landsat 8(USGS, 2023a, ) data to achieve one LSTreconstruction every 16 days.Recently, the Landsat program launched a companion satellite, Landsat9(USGS, 2023b, ), carrying a nearly identical TIRS as Landsat 8.Landsat 8 and Landsat 9 are phased eight days apart.By incorporating data from both satellites, we can reduce the time gap between consecutive satellite visits to just eight days.This enables us to capture LST measurements at a higher temporal resolution.

Another avenue to increasing the temporal resolution is to use satelliteproducts with shorter revisit cycles. For example, theSentinel-3(Donlon etal.,, 2012) program provides a temporal resolution ofat least once per day (at the equator). Unfortunately, the spatial resolutionof their LST products is lower than that of Landsat 8 (e.g.,Sentinel-3 at$1\text{\,}\mathrm{k}\mathrm{m}$).Additional benchmarking on the selected data source is required before applying ISLAND, as performance could vary based on spatio-temporal resolution.

Incorporating deep learning:The basis of our model is a set of filters designed around adjacency and temporal properties of thermal signatures.As seen in Algorithms 1 and 2, these filters are hand crafted to explicitly represent these relationships.While we clearly demonstrated the effectiveness of ISLAND through qualitative and quantitative analysis, we acknowledge that a well-designed deep learning algorithm has the potential to achieve even better performance.Here, we highlight a few examples.Firstly, deep learning models have the capability to capture complex inter-class relationships between different land cover labels.This could enhance the overall accuracy of our interpolator.Secondly, a dynamic spatial channel that adapts based on occlusion characteristics could be incorporated, allowing the model to better handle varying cloud cover conditions.Additionally, an optimized weighting scheme, an improved cloud detection filter, and an updated NLCD land cover dataset to account for changes in land cover could be integrated into a deep learning framework.Lastly, integrating additional data sources, such as other satellite data andground-based observations, could further reduce reconstruction errors, andthere are existing deep learning techniques(Han etal.,, 2024) to perform datafusion.Given the complexity of the problem and the non-linear nature ofLST(Wu etal.,, 2021), deep learning is a suitable direction for futurework, but designing and training a deep learning framework might requireextensive research.

6.4 Advancing Existing State-of-the-Art LST Estimates

In this paper, we showed that a large fraction of Landsat measurements are occluded by clouds. As a consequence, the actual usable temporal resolution of Landsat is significantly reduced, falling below once per month when subject to frequent cloud occlusions. The role of ISLAND is to mitigate cloud contamination in LST images, maximizing its usable temporal resolution.

Indeed, the addition of ISLAND represents an advance over existing LST products, as shown in Fig.10.Generally speaking, there is a trade-off between spatial resolution and temporal resolution for satellite LST products.Amongst all available satellite measurements, Landsat 8 (along with the later-launched Landsat 9) offers the best spatial resolution at $30\text{\,}\mathrm{m}$, with 16-day revisit cycles.Operating in sun-synchronous orbits (SSO), Landsat 8 has the advantage of providing global coverage and maintaining time-constant illumination conditions of the observed surfaces (except for seasonal variations).The MODIS program(Wan,, 2013), consisting of a pair of satellites namedAqua and Terra, offers a much higher resolution at daily revisit cycles butprovides LST data at a much lower spatial resolution of $1\text{\,}\mathrm{km}$.Satellites in geostationary orbits do provide higher temporal resolution at the expense of spatial resolution and global coverage.For example, GridSat-B1(Knapp etal.,, 2011; Knapp,, 2014) provides brightnesstemperature (BT) data at a resolution of $7792\text{\,}\mathrm{m}$.

In addition to satellite-based methods, climate reanalysis data provides an alternative approach to obtaining LST.These products are generally designed to maintain the best possible physicaland temporal consistency and require prohibitive computationalcosts(Hakim etal.,, 2016).The spatial resolution of these products is not comparable to satellite-basedmethods; HRRR(James etal.,, 2022) provides climate data at $3\text{\,}\mathrm{km}$ andERA5(MuñozSabater,, 2019) at around $11\text{\,}\mathrm{km}$.Fig.10(b) provides a visual comparison of the spatial resolution of LST from ERA5 and ISLAND.Due to relatively low spatial resolution, the urban spatial structure over Philadelphia is indistinguishable in ERA5 skin temperature fields. In contrast, our method effectively removes cloud contamination and produces a high-resolution reconstruction of LST at $30\text{\,}\mathrm{m}$ resolution.

As shown in Fig.5.5, many downstream applications generally benefit from increased spatial and temporal resolution.For dense urban settings, high spatial resolution is particularly desirable, making the Landsat data a preferred choice.In Fig.6.3, we discussed the potential of incorporating measurements from Landsat 9 to further enhance the temporal resolution to 8 days.This advancement would bring us closer to achieving consistent weekly measurements at a spatial resolution of $30\text{\,}\mathrm{m}$.Such a combination of high spatial and temporal resolution of LST data is instrumental to our understanding of urban areas.

6.5 Conclusions

This paper introduces ISLAND, a novel model designed to address the issue of cloud occlusion in satellite LST images. ISLAND removes occlusion by estimating LST pixel values through a set of spatio-temporal filters. These filters account for the land cover class, resulting in higher LST reconstruction accuracy. ISLAND addresses a fundamental limitation of LST retrieval via remote sensing, thereby dramatically increasing the number of serviceable LST images via a robust mechanism to mitigate cloud contamination. These improvements enable nearly bi-weekly coverage of LST at $30\text{\,}\mathrm{m}$ resolution over the CONUS region, a large advance over previously available LST products derived from remote sensing. We show ISLAND can operate in a variety of land cover types and cloud occlusion scenarios in both simulations and in situ evaluations. Overall, ISLAND provides a promising framework for a multitude of scientific applications that require high-resolution, frequent observations of LST, including but not limited to (1) urban heat island effects, (2) derivation of surface temperature trends, and (3) social vulnerability and urban heat stress.