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

Yuhao Liu
Rice University
yuhao.liu@rice.edu&Pranavesh Panakkal
Rice University
pranavesh@rice.edu&Sylvia Dee
Rice University
sylvia.dee@rice.edu&Guha Balakrishnan
Rice University
guha@rice.edu&Jamie Padgett
Rice University
jamie.padgett@rice.edu&Ashok Veeraraghavan
Rice University
vashok@rice.edu
Corresponding author

Abstract

Cloud occlusion is a common problem in the field of remote sensing, particularly for retrieving Land Surface Temperature (LST). Remote sensing thermal instruments onboard operational satellites are supposed to enable frequent and high-resolution observations over land; unfortunately, clouds adversely affect thermal signals by blocking outgoing longwave radiation emission from the Earth’s surface, interfering with the retrieved ground emission temperature. Such cloud contamination severely reduces the set of serviceable LST images for downstream applications, making it impractical to perform intricate time-series analysis of LST. In this paper, we introduce a novel method to remove cloud occlusions from Landsat 8 LST images. We call our method ISLAND, an acronym for Interpolating Land Surface Temperature using land cover. Our approach uses LST images from Landsat 8 (at 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG resolution with 16-day revisit cycles) and the NLCD land cover dataset. Inspired by Tobler’s first law of Geography, ISLAND predicts occluded LST through a set of spatio-temporal filters that perform distance-weighted spatio-temporal interpolation. A critical feature of ISLAND is that the filters are land cover-class aware, making it particularly advantageous in complex urban settings with heterogeneous land cover types and distributions. Through qualitative and quantitative analysis, we show that ISLAND achieves robust reconstruction performance across a variety of cloud occlusion and surface land cover conditions, and with a high spatio-temporal resolution. We provide a public dataset of 20 U.S. cities with pre-computed ISLAND LST outputs. Using several case studies, we demonstrate that ISLAND opens the door to a multitude of high-impact urban and environmental applications across the continental United States.

Keywords cloud removal \cdotland surface temperature \cdotthermal imaging \cdotLandsat \cdotland cover

1 Introduction

Land surface temperature (LST) is a fundamental aspect of Earth’s climate system, it attracts extensive studies across diverse disciplines. The wide array of fields include climate change (Horton etal.,, 2016; Seneviratne etal.,, 2021), urban planning (Sobrino etal.,, 2013; Huang and Wang,, 2019; Osborne and Alvares-Sanches,, 2019), vegetation and land cover changes(Gomez-Martinez etal.,, 2021; Chun and Guldmann,, 2018), and human health (Orimoloye etal.,, 2018; Dee etal.,, 2022; Dugord etal.,, 2014). LST changes rapidly both in space and time due to the strong heterogeneity of land surface characteristics and the short timescales of weather. As a consequence, accurate characterization of LST requires dense spatial and temporal sampling (Li etal.,, 2013).

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

Growing recognition of the importance of accurate LST observations has drivenrapid advances in remote sensing(Li etal.,, 2013).Satellite-based thermal infrared (TIR) data provides LST measurements with high spatial and temporal resolution at a global scale.As an example, the National Aeronautics and Space Administration (NASA) Landsat8 satellite was launched in 2013, and its platform has provided detailed LSTdata over the last decade(USGS, 2023a, ; USGS, 2023b, ).The Landsat 8 and 9 satellites each include a Thermal Infrared Sensor (TIRS) onboard, which collects LST data at a spatial resolution of 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG with a revisit cycle of 16 days.The high spatial resolution of TIRS has enabled numerous studies of LST overheterogeneous regions like urbancenters(Huang and Wang,, 2019; Sobrino etal.,, 2013; Streutker,, 2003)and wetlands(Eisavi etal.,, 2016; Demarquet etal.,, 2023).

Despite these advancements, cloud occlusion persists as a substantial obstacle to achieving reliable spaceborne retrievals of LST.Clouds adversely affect LST readings by blocking thermal radiation emitted from Earth’s surface.This results in cloud-contaminated pixels exhibiting considerably lower valuescompared to their true values, rendering the affected imagesunusable(Jin etal.,, 2019).Unfortunately, cloud contamination is a frequent phenomenon in Landsat images,with studies indicating that an average of 35% of Landsat images globallycontain missing data due to cloudcontamination(Roy etal.,, 2008).Moreover, most cities experience enhanced daytime cloud cover compared to theirsurrounding rural regions(Vo etal.,, 2023), making cloud-free LST imageseven harder to obtain for cities.Furthermore, regions at higher latitudes, such as Pacific Northwest cities likePortland and Seattle, are prone to rainier weather conditions,leading to more frequent cloud contaminations(NOAA,, 2020).Under these circ*mstances, obtaining a single cloud-free image may be infeasible for weeks or even months, severely limiting the utility of Landsat LST data.Collectively, these barriers have limited practical high spatial and temporal sampling of LST. Thus, the temporal dynamics of LST over complex, heterogeneous terrains remain under-observed, under-studied, and under-constrained, especially at scale.

In this paper, we present a novel method to mitigate the effects of cloud contamination in satellite LST images.Our method incorporates Tobler’s First Law of Geography (TFL), which statesthat “everything is related to everything else, but near things are morerelated than distant things”(Tobler,, 1970).In addition to the distance-decay effect from TFL, we also incorporate multi-temporal information and, crucially, integrate land cover data into our model.Land cover data contains information about physical land types, such as forests, open water, and urban/developed areas.While existing studies have demonstrated the strong relationship between landcover types andLST(Chaudhuri and Mishra,, 2016; Zhao etal.,, 2020; Imran etal.,, 2021), ourmodel represents the first attempt to utilize land cover to infer occluded LSTpixel values in satellite images.We call our model ISLAND, an acronym for Interpolating Land Surface Temperature using land cover.ISLAND performs interpolation using a set of spatio-temporal filters to capture surrounding pixel values and historical patterns.Notably, these filters are designed to be sensitive to land cover classes, enabling class-specific pattern capture and higher reconstruction accuracy.

We demonstrate that ISLAND provides robust LST reconstruction performance across various occlusion and land cover conditions.Our results indicate that ISLAND greatly improves the practical temporal resolution of Landsat LST images by robustly estimating cloud-contaminated LST pixels.Simulation evaluation shows that the RMSE in our reconstructed Landsat 8 LST images is around 1.65limit-from1.651.65-1.65 -2.62Ktimes2.62kelvin2.62\text{\,}\mathrm{K}start_ARG 2.62 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG. In situ evaluation shows that our RMSE is approximately 3.90limit-from3.903.90-3.90 -4.75Ktimes4.75kelvin4.75\text{\,}\mathrm{K}start_ARG 4.75 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG.We present three illustrative examples underscoring the utility of ISLAND, namely, (1) urban heat island effects, (2) derivation of surface temperature trends, and (3) social vulnerability and urban heat stress.

ISLAND leverages publicly available data from Landsat 8(USGS, 2023a, ) andNational Land Cover Database (NLCD)(Yang etal.,, 2018). This approach ensures the accessibility andtransferability of ISLAND, as it is open-source and can be easily deployed inany region within the continental United States (CONUS) as per userrequirements.All source code111https://github.com/Way-Yuhao/ISLAND is available to the public, accompanied by a dataset222https://doi.org/10.17603/ds2-3rf5-sd58 comprising 20 urban regions with pre-computed ISLAND LST outputs.We envision that ISLAND will provide tremendous operational value for a variety of applications in Earth, geospatial, and social sciences, and ultimately pave the way for a new generation of LST studies using remote sensing.

2 Related Work

The refinement of cloud-removing algorithms is an open area of research in the field of remote sensing.Several existing studies develop and discuss cloud removal algorithms for remotely acquired LST.Unfortunately, these algorithms generally target lower spatial resolutions andhom*ogeneous land cover types(Wu etal.,, 2021).For example,Yu etal., (2019) proposed a method to reconstruct ModerateResolution Imaging Spectroradiometer (MODIS) LST over cloudy pixels using landenergy balance theory and similar pixels at a spatial resolution of1kmtimes1kilometer1\text{\,}\mathrm{km}start_ARG 1 end_ARG start_ARG times end_ARG start_ARG roman_km end_ARG.Zeng etal., (2015) proposed a spatiotemporal technique to reconstruct MODIS LSTproducts using regression and multispectral ancillary data to classify pixels,which improves reconstruction accuracy.Unfortunately, MODIS LST at 1kmtimes1kilometer1\text{\,}\mathrm{km}start_ARG 1 end_ARG start_ARG times end_ARG start_ARG roman_km end_ARG is too coarse to resolve urban infrastructures, such as buildings and roads.

There are existing studies that reconstruct cloudy-sky LST at the spatial resolution of Landsat 8 pixels (30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG).Wang etal., (2019) proposed a method to reconstruct cloud-sky Landsat 8LST by considering solar-cloud-satellite geometry.However, this method requires a temporally adjacent clear-sky image as a reference, which is hard to obtain, limiting its robustness.Furthermore, the accuracy of their reconstructed Landsat LST was not validated against in situ LST measurements.Zhu etal., (2022) reconstruct Landsat 8 LST data using an annualtemperature cycle (ATC) model and adjacent spatial information from similarpixels.Zhu etal., (2022) validated their results against in situ LST measurementsat six Surface Radiation Budget Network (SURFRAD) sites(Augustine etal., 2000a, ).However, there are several limitations to their study.First, their reconstruction method may fail when the cloud cover percentage exceeds 70%, limiting the scope of their application.Additionally, the ATC model assumes a constant mean annual surface temperature per region, suggesting that their model likely underestimates LST under global warming and the increased prevalence of the urban heat island effect.Lastly, the primary limitation forZhu etal., (2022) and all existingstudies mentioned in this section is that these studies are mainly designed andvalidated on relatively hom*ogeneous regions of land cover, such as cropland,shrubland, and grassland; their results were not demonstrated over urbanregions.This study addresses this gap and provides a scalable method to generate high-resolution LST, even for highly urbanized regions with complex temperature feedback and heterogeneous urban land surface types.We achieve this by explicitly employing NLCD land cover labels and using the satellite LST product with the highest resolution (Landsat 8) to maximize its usability in cities.

Our model shares some of the high-level design principles with existing models, such as the use of multitemporal remote sensing data.To the best of our knowledge, ISLAND is the first algorithm that incorporates land cover labels as ancillary data for cloud removal in LST images.We summarize our key contributions as follows:(1) We use NLCD land cover labels to accurately reconstruct Landsat 8 LST under cloudy-sky conditions.(2) Using spatial adjacency and multi-temporal filters, ISLAND effectively removes cloud contaminations from Landsat 8 LST images even under severe occlusion, thereby dramatically improving the temporal resolution of Landsat 8 LST products.(3) Simulation evaluation shows that the RMSE in our reconstructed Landsat 8 LST images is around 1.65limit-from1.651.65-1.65 -2.62Ktimes2.62kelvin2.62\text{\,}\mathrm{K}start_ARG 2.62 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG. In situ evaluation shows that our RMSE is approximately 3.90limit-from3.903.90-3.90 -4.75Ktimes4.75kelvin4.75\text{\,}\mathrm{K}start_ARG 4.75 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG.(4) We release a public dataset of cloud-free, reconstructed LST maps for 20 US cities from 2019–2023.

3 Data

In this section, we explain the implementation of the data compilation steps of our model. See Fig.1 for a visual overview.

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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG 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 approximately2.2Ktimes2.2kelvin2.2\text{\,}\mathrm{K}start_ARG 2.2 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG(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 T~H×W~𝑇superscript𝐻𝑊\tilde{T}\in\mathbb{R}^{H\times W}over~ start_ARG italic_T end_ARG ∈ blackboard_R start_POSTSUPERSCRIPT italic_H × italic_W end_POSTSUPERSCRIPT, where H𝐻Hitalic_H and W𝑊Witalic_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{0,1}H×W𝑂superscript01𝐻𝑊O\in\{0,1\}^{H\times W}italic_O ∈ { 0 , 1 } start_POSTSUPERSCRIPT italic_H × italic_W end_POSTSUPERSCRIPT 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 θ[0,1]𝜃01\theta\in[0,1]italic_θ ∈ [ 0 , 1 ] that measures the fraction of pixels occluded.Formally, we have θ=𝒑O𝒑HW𝜃subscript𝒑subscript𝑂𝒑𝐻𝑊\theta=\frac{\sum_{\boldsymbol{p}}O_{\boldsymbol{p}}}{H\cdot W}italic_θ = divide start_ARG ∑ start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT italic_O start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT end_ARG start_ARG italic_H ⋅ italic_W end_ARG, where Opsubscript𝑂𝑝O_{p}italic_O start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT denotes a pixel of O𝑂Oitalic_O.θ𝜃\thetaitalic_θ is an important metric that measures the severity of cloud contamination.A higher θ𝜃\thetaitalic_θ 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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG resolution with a release cycle of once per 3years.We denote the NLCD land cover image as LH×W𝐿superscript𝐻𝑊L\in\mathbb{Z}^{H\times W}italic_L ∈ blackboard_Z start_POSTSUPERSCRIPT italic_H × italic_W end_POSTSUPERSCRIPT.

3.3 In situ data

We collect in situ LST data at four Surface Radiation Budget Network (SURFRAD)sites(Augustine etal., 2000b, ).We use SURFRAD in situ data to validate our reconstruction results under cloudy conditions.The RMSE of SURFRAD in situ LST measurements is approximately0.5–0.8Ktimes0.8kelvin0.8\text{\,}\mathrm{K}start_ARG 0.8 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG(Wang and Liang,, 2009).SURFRAD measures downwelling and upwelling longwave flux every minute.The Stefan–Boltzmann law states

F=(1ϵ)F+ϵbσTs4superscript𝐹1italic-ϵsuperscript𝐹subscriptitalic-ϵ𝑏𝜎superscriptsubscript𝑇𝑠4F^{\uparrow}=(1-\epsilon)F^{\downarrow}+\epsilon_{b}\sigma T_{s}^{4}italic_F start_POSTSUPERSCRIPT ↑ end_POSTSUPERSCRIPT = ( 1 - italic_ϵ ) italic_F start_POSTSUPERSCRIPT ↓ end_POSTSUPERSCRIPT + italic_ϵ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT italic_σ italic_T start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT(1)

where Fsuperscript𝐹F^{\uparrow}italic_F start_POSTSUPERSCRIPT ↑ end_POSTSUPERSCRIPT and Fsuperscript𝐹F^{\downarrow}italic_F start_POSTSUPERSCRIPT ↓ end_POSTSUPERSCRIPT are measured upwelling and downwelling longwave flux, respectively, ϵbsubscriptitalic-ϵ𝑏\epsilon_{b}italic_ϵ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT is the broadband longwave surface emissivity, σ𝜎\sigmaitalic_σ is the Stefan–Boltzmann constant, and Tssubscript𝑇𝑠T_{s}italic_T start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT is the surface skin temperature (equivalent to LST).

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

Ts={1εbσ[F(1ϵ)F]}0.25subscript𝑇𝑠superscript1subscript𝜀𝑏𝜎delimited-[]superscript𝐹1italic-ϵsuperscript𝐹0.25T_{s}=\left\{\frac{1}{\varepsilon_{b}\sigma}\left[F^{\uparrow}-(1-\epsilon)F^{%\downarrow}\right]\right\}^{0.25}italic_T start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = { divide start_ARG 1 end_ARG start_ARG italic_ε start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT italic_σ end_ARG [ italic_F start_POSTSUPERSCRIPT ↑ end_POSTSUPERSCRIPT - ( 1 - italic_ϵ ) italic_F start_POSTSUPERSCRIPT ↓ end_POSTSUPERSCRIPT ] } start_POSTSUPERSCRIPT 0.25 end_POSTSUPERSCRIPT(2)

FollowingMalakar etal., (2018), we estimate the broadband emissivity(ϵbsubscriptitalic-ϵ𝑏\epsilon_{b}italic_ϵ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT) via a spectral-to-broadband regressionrelationship(Ogawa etal.,, 2008) using emissivity values from ASTERGEDv3 product(Hulley etal.,, 2015)

ϵb=0.128+0.014ϵa10+0.145ϵa11+0.241ϵa12+0.467ϵa13+0.004ϵa14subscriptitalic-ϵ𝑏0.1280.014superscriptsubscriptitalic-ϵ𝑎100.145superscriptsubscriptitalic-ϵ𝑎110.241superscriptsubscriptitalic-ϵ𝑎120.467superscriptsubscriptitalic-ϵ𝑎130.004superscriptsubscriptitalic-ϵ𝑎14\epsilon_{b}=0.128+0.014\epsilon_{a}^{10}+0.145\epsilon_{a}^{11}+0.241\epsilon%_{a}^{12}+0.467\epsilon_{a}^{13}+0.004\epsilon_{a}^{14}italic_ϵ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT = 0.128 + 0.014 italic_ϵ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 10 end_POSTSUPERSCRIPT + 0.145 italic_ϵ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 11 end_POSTSUPERSCRIPT + 0.241 italic_ϵ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT + 0.467 italic_ϵ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 13 end_POSTSUPERSCRIPT + 0.004 italic_ϵ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 14 end_POSTSUPERSCRIPT(3)

where ϵajsuperscriptsubscriptitalic-ϵ𝑎𝑗\epsilon_{a}^{j}italic_ϵ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_j end_POSTSUPERSCRIPT, j=1014𝑗1014j=10-14italic_j = 10 - 14 are narrowband ASTER emissivities centered on 8.3, 8.6, 9.1, 10.6, and 11.3µmtimes11.3micrometer11.3\text{\,}\mathrm{\SIUnitSymbolMicro m}start_ARG 11.3 end_ARG start_ARG times end_ARG start_ARG roman_µ roman_m end_ARG, respectively.

4 Methods

After acquiring the required inputs from Sec.3, we perform interpolation to predict cloud-occluded pixels.Our interpolator uses two complementary mechanisms: a spatial channel and a temporal channel.

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

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𝐿Litalic_L and occlusion O𝑂Oitalic_O).

Consider our goal to be predicting T~~𝑇\tilde{T}over~ start_ARG italic_T end_ARG at an occluded pixel location 𝒑𝒑\boldsymbol{p}bold_italic_p.Under a low occlusion factor θ𝜃\thetaitalic_θ, we estimate T~𝒑subscript~𝑇𝒑\tilde{T}_{\boldsymbol{p}}over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT 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.

We use a 2D Gaussian filter to model the distance-decay effect.Consider a pixel 𝒒𝒒\boldsymbol{q}bold_italic_q within a f×f𝑓𝑓f\times fitalic_f × italic_f neighborhood N𝑁Nitalic_N centered at 𝒑𝒑\boldsymbol{p}bold_italic_p.Let x=𝒑𝒒𝑥delimited-∥∥𝒑𝒒x=\lVert\boldsymbol{p}-\boldsymbol{q}\rVertitalic_x = ∥ bold_italic_p - bold_italic_q ∥ be the Euclidean distance between 𝒑𝒑\boldsymbol{p}bold_italic_p and 𝒒𝒒\boldsymbol{q}bold_italic_q.We compute the proximity weight Gσ(x)subscript𝐺𝜎𝑥G_{\sigma}(x)italic_G start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ( italic_x ) using a Gaussian kernel with standard deviation σ=f/2𝜎𝑓2\sigma=f/2italic_σ = italic_f / 2:

Gσ(x)=12πσ2exp(x22σ2)subscript𝐺𝜎𝑥12𝜋superscript𝜎2expsuperscript𝑥22superscript𝜎2G_{\sigma}(x)=\frac{1}{2\pi\sigma^{2}}\text{exp}\biggl{(}-\frac{x^{2}}{2\sigma%^{2}}\biggr{)}italic_G start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ( italic_x ) = divide start_ARG 1 end_ARG start_ARG 2 italic_π italic_σ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG exp ( - divide start_ARG italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_σ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG )(4)

Gσ(x)subscript𝐺𝜎𝑥G_{\sigma}(x)italic_G start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ( italic_x ) decreases the influence of distance pixels while prioritizing the influence of nearby pixels.

We further modify the Gaussian filtering so that we only consider 𝒒𝒒\boldsymbol{q}bold_italic_q that is cloud-free and has the same land cover class as 𝒑𝒑\boldsymbol{p}bold_italic_p.Let O¯¯𝑂\bar{O}over¯ start_ARG italic_O end_ARG be the inverse of the occlusion mask O𝑂Oitalic_O, so that O¯𝒒=1subscript¯𝑂𝒒1\bar{O}_{\boldsymbol{q}}=1over¯ start_ARG italic_O end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT = 1 if there is no cloud contamination at 𝒒𝒒\boldsymbol{q}bold_italic_q, and O¯𝒒=0subscript¯𝑂𝒒0\bar{O}_{\boldsymbol{q}}=0over¯ start_ARG italic_O end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT = 0 otherwise.To constrain land cover class, let function fLsubscript𝑓𝐿f_{L}italic_f start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT evaluate the land cover class of two pixels 𝒑𝒑\boldsymbol{p}bold_italic_p and 𝒒𝒒\boldsymbol{q}bold_italic_q such that fL(𝒑,𝒒)=1subscript𝑓𝐿𝒑𝒒1f_{L}(\boldsymbol{p},\boldsymbol{q})=1italic_f start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( bold_italic_p , bold_italic_q ) = 1 if the corresponding land cover classes of 𝒑𝒑\boldsymbol{p}bold_italic_p and 𝒒𝒒\boldsymbol{q}bold_italic_q are the same (i.e.,L𝒑=L𝒒subscript𝐿𝒑subscript𝐿𝒒L_{\boldsymbol{p}}=L_{\boldsymbol{q}}italic_L start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT = italic_L start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT), and fL(𝒑,𝒒)=0subscript𝑓𝐿𝒑𝒒0f_{L}(\boldsymbol{p},\boldsymbol{q})=0italic_f start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( bold_italic_p , bold_italic_q ) = 0 otherwise.Formally, we write our weighted average local filter as follows:

T~𝒑=1α𝒒NGσ(𝒑𝒒)O¯𝒒fL(𝒑,𝒒)T~𝒒subscript~𝑇𝒑1𝛼subscript𝒒𝑁subscript𝐺𝜎delimited-∥∥𝒑𝒒subscript¯𝑂𝒒subscript𝑓𝐿𝒑𝒒subscript~𝑇𝒒\tilde{T}_{\boldsymbol{p}}=\frac{1}{\alpha}\sum_{\boldsymbol{q}\in N}G_{\sigma%}(\lVert\boldsymbol{p}-\boldsymbol{q}\rVert)~{}\bar{O}_{\boldsymbol{q}}~{}f_{L%}(\boldsymbol{p},\boldsymbol{q})~{}\tilde{T}_{\boldsymbol{q}}over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG italic_α end_ARG ∑ start_POSTSUBSCRIPT bold_italic_q ∈ italic_N end_POSTSUBSCRIPT italic_G start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ( ∥ bold_italic_p - bold_italic_q ∥ ) over¯ start_ARG italic_O end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( bold_italic_p , bold_italic_q ) over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT(5)

where α𝛼\alphaitalic_α is a normalization parameter that ensures weights sum to 1111 within each neighborhood (i.e., 𝒒Nfor-all𝒒𝑁\forall\boldsymbol{q}\in N∀ bold_italic_q ∈ italic_N):

α=𝒒NGσ(𝒑𝒒)O¯𝒒fL(𝒑,𝒒)𝛼subscript𝒒𝑁subscript𝐺𝜎delimited-∥∥𝒑𝒒subscript¯𝑂𝒒subscript𝑓𝐿𝒑𝒒\alpha=\sum_{\boldsymbol{q}\in N}G_{\sigma}(\lVert\boldsymbol{p}-\boldsymbol{q%}\rVert)~{}\bar{O}_{\boldsymbol{q}}~{}f_{L}(\boldsymbol{p},\boldsymbol{q})italic_α = ∑ start_POSTSUBSCRIPT bold_italic_q ∈ italic_N end_POSTSUBSCRIPT italic_G start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ( ∥ bold_italic_p - bold_italic_q ∥ ) over¯ start_ARG italic_O end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( bold_italic_p , bold_italic_q )(6)

Local filtering (defined in Eq.(5)) works well when the occlusion factor is low.For a high θ𝜃\thetaitalic_θ, there are fewer neighboring pixels available, and local filtering leads to noisy or even invalid estimations.Instead of local filtering, we resort to global averaging when encountering high θ𝜃\thetaitalic_θ.Here we estimate the pixel value of 𝒑𝒑\boldsymbol{p}bold_italic_p with the average temperature of all non-occluded pixels that have the same land cover class, that is

μc=1T𝒒T~{𝒑}O¯𝒒fL(𝒑,𝒒)T~𝒒subscript𝜇𝑐1𝑇subscript𝒒~𝑇𝒑subscript¯𝑂𝒒subscript𝑓𝐿𝒑𝒒subscript~𝑇𝒒\mu_{c}=\frac{1}{T}\sum_{\boldsymbol{q}\in\tilde{T}\setminus\{\boldsymbol{p}\}%}\bar{O}_{\boldsymbol{q}}~{}f_{L}(\boldsymbol{p},\boldsymbol{q})~{}\tilde{T}_{%\boldsymbol{q}}italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG italic_T end_ARG ∑ start_POSTSUBSCRIPT bold_italic_q ∈ over~ start_ARG italic_T end_ARG ∖ { bold_italic_p } end_POSTSUBSCRIPT over¯ start_ARG italic_O end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( bold_italic_p , bold_italic_q ) over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT(7)

where c𝑐citalic_c is the land cover class of the occluded pixel 𝒑𝒑\boldsymbol{p}bold_italic_p, and T𝑇Titalic_T is the total number of non-occluded pixels with the same land cover class as 𝒑𝒑\boldsymbol{p}bold_italic_p, excluding 𝒑𝒑\boldsymbol{p}bold_italic_p itself.Contrary to local filtering (Eq.(5)), here the averaging is performed across the entire image, rather than across some local f×f𝑓𝑓f\times fitalic_f × italic_f window.As a result, we are no longer able to capture proximity effects.Therefore, outputs tend to be blurry due to spatial averaging.

Algorithm 1 shows the implementation of the spatial channel.Let T~~𝑇\tilde{T}over~ start_ARG italic_T end_ARG be the occluded LST image as input, f𝑓fitalic_f be the size of the local window, and θsuperscript𝜃\theta^{*}italic_θ start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT be the maximum occlusion factor threshold for local filtering.We obtain interpolated image T^spsubscript^𝑇𝑠𝑝\hat{T}_{sp}over^ start_ARG italic_T end_ARG start_POSTSUBSCRIPT italic_s italic_p end_POSTSUBSCRIPT as follows:

1:T~,f,θ~𝑇𝑓superscript𝜃\tilde{T},f,\theta^{*}over~ start_ARG italic_T end_ARG , italic_f , italic_θ start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT

2:ifθ<θ𝜃superscript𝜃\theta<\theta^{*}italic_θ < italic_θ start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPTthen \triangleright local filtering

3:forT~𝒑T~subscript~𝑇𝒑~𝑇\tilde{T}_{\boldsymbol{p}}\in\tilde{T}over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT ∈ over~ start_ARG italic_T end_ARGdo

4:ifO𝒑=1subscript𝑂𝒑1O_{\boldsymbol{p}}=1italic_O start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT = 1then

5:T^𝒑1α𝒒NGσ(𝒑𝒒)O¯𝒒fL(𝒑,𝒒)T~𝒒subscript^𝑇𝒑1𝛼subscript𝒒𝑁subscript𝐺𝜎delimited-∥∥𝒑𝒒subscript¯𝑂𝒒subscript𝑓𝐿𝒑𝒒subscript~𝑇𝒒\hat{T}_{\boldsymbol{p}}\leftarrow\frac{1}{\alpha}\sum_{\boldsymbol{q}\in N}G_%{\sigma}(\lVert\boldsymbol{p}-\boldsymbol{q}\rVert)~{}\bar{O}_{\boldsymbol{q}}%~{}f_{L}(\boldsymbol{p},\boldsymbol{q})~{}\tilde{T}_{\boldsymbol{q}}over^ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT ← divide start_ARG 1 end_ARG start_ARG italic_α end_ARG ∑ start_POSTSUBSCRIPT bold_italic_q ∈ italic_N end_POSTSUBSCRIPT italic_G start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ( ∥ bold_italic_p - bold_italic_q ∥ ) over¯ start_ARG italic_O end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( bold_italic_p , bold_italic_q ) over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_q end_POSTSUBSCRIPT

6:endif

7:endfor

8:else\triangleright global averaging

9:forT~𝒑T~subscript~𝑇𝒑~𝑇\tilde{T}_{\boldsymbol{p}}\in\tilde{T}over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT ∈ over~ start_ARG italic_T end_ARGdo

10:ifO𝒑=1subscript𝑂𝒑1O_{\boldsymbol{p}}=1italic_O start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT = 1then

11:cL𝒑𝑐subscript𝐿𝒑c\leftarrow L_{\boldsymbol{p}}italic_c ← italic_L start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT

12:T^𝒑μcsubscript^𝑇𝒑subscript𝜇𝑐\hat{T}_{\boldsymbol{p}}\leftarrow\mu_{c}over^ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT ← italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT

13:endif

14:endfor

15:endif

16:return T^spT^superscript^𝑇𝑠𝑝^𝑇\hat{T}^{sp}\leftarrow\hat{T}over^ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT italic_s italic_p end_POSTSUPERSCRIPT ← over^ start_ARG italic_T end_ARG

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

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

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 T~~𝑇\tilde{T}over~ start_ARG italic_T end_ARG and other available images T~superscript~𝑇\tilde{T}^{\prime}over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT.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 δ𝛿\deltaitalic_δ 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 θ𝜃\thetaitalic_θ.We have observed that minimally occluded reference frames lead to smaller interpolation errors.Therefore we only consider selecting reference frames whose θ𝜃\thetaitalic_θ is below a certain maximum tolerable threshold, θmaxsubscript𝜃max\theta_{\text{max}}italic_θ start_POSTSUBSCRIPT max end_POSTSUBSCRIPT (see red dotted line in Fig.3(b)).We call the set of images satisfying these two conditions candidate reference frames 𝒯refsubscriptsuperscript𝒯𝑟𝑒𝑓\mathcal{T}^{\ast}_{ref}caligraphic_T start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT.Within 𝒯refsubscriptsuperscript𝒯𝑟𝑒𝑓\mathcal{T}^{\ast}_{ref}caligraphic_T start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT, we prioritize selecting images that are captured closer to the target frame in time, as measured by the temporal difference Δ(|T~T~|)Δsuperscript~𝑇~𝑇\Delta(|\tilde{T}^{\prime}-\tilde{T}|)roman_Δ ( | over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT - over~ start_ARG italic_T end_ARG | ).To achieve this, we choose a subset of n𝑛nitalic_n frames from 𝒯refsubscriptsuperscript𝒯𝑟𝑒𝑓\mathcal{T}^{\ast}_{ref}caligraphic_T start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT with the lowest Δ(|T~T~|)Δsuperscript~𝑇~𝑇\Delta(|\tilde{T}^{\prime}-\tilde{T}|)roman_Δ ( | over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT - over~ start_ARG italic_T end_ARG | ).We call this selected subset of reference frames 𝒯refsubscript𝒯𝑟𝑒𝑓\mathcal{T}_{ref}caligraphic_T start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT.

Spatial channel pre-processing:Given a selection of reference frames 𝒯refsubscript𝒯𝑟𝑒𝑓\mathcal{T}_{ref}caligraphic_T start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT, 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 θmaxsubscript𝜃𝑚𝑎𝑥\theta_{max}italic_θ start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT means that the spatial channel only needs to interpolate a minimal amount of occlusion.

Linear adjustments:After pre-processing 𝒯refsubscript𝒯𝑟𝑒𝑓\mathcal{T}_{ref}caligraphic_T start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT, we apply a linear adjustment to all pixels in each class. Fig.3(a) shows that the changes in the mean LST (denoted as μcsubscript𝜇𝑐\mu_{c}italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT) 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 μcsubscript𝜇𝑐\mu_{c}italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT (i.e.,ΔμcΔsubscript𝜇𝑐\Delta\mu_{c}roman_Δ italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT) between two dates to all pixels belonging to the corresponding land cover class.ΔμcΔsubscript𝜇𝑐\Delta\mu_{c}roman_Δ italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT visually translates to the length of each bar in Fig.3(a).We repeat this process for all images in 𝒯refsubscript𝒯𝑟𝑒𝑓\mathcal{T}_{ref}caligraphic_T start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT.

Reference frame-based interpolation:After the previous two steps, we obtain a set of linearly-adjusted reference frames.The interpolated LST image, T^tempsuperscript^𝑇𝑡𝑒𝑚𝑝\hat{T}^{temp}over^ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT italic_t italic_e italic_m italic_p end_POSTSUPERSCRIPT, 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𝑛nitalic_n be the number of reference frames to be selected, δ𝛿\deltaitalic_δ be the temporal bracket duration, and θmaxsubscript𝜃𝑚𝑎𝑥\theta_{max}italic_θ start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT be the maximum tolerable occlusion factor. We compute the interpolated image T^tempsubscript^𝑇𝑡𝑒𝑚𝑝\hat{T}_{temp}over^ start_ARG italic_T end_ARG start_POSTSUBSCRIPT italic_t italic_e italic_m italic_p end_POSTSUBSCRIPT as follows:

1:T~,n,δ,θmax~𝑇𝑛𝛿subscript𝜃𝑚𝑎𝑥\tilde{T},n,\delta,\theta_{max}over~ start_ARG italic_T end_ARG , italic_n , italic_δ , italic_θ start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT

2:𝒯ref={T~𝒯|δT~<δθT~<θmax}subscriptsuperscript𝒯𝑟𝑒𝑓conditional-setsuperscript~𝑇𝒯subscript𝛿superscript~𝑇𝛿subscript𝜃superscript~𝑇subscript𝜃𝑚𝑎𝑥\mathcal{T}^{\ast}_{ref}=\{\tilde{T}^{\prime}\in\mathcal{T}|\delta_{\tilde{T}^%{\prime}}<\delta\land\theta_{\tilde{T}^{\prime}}<\theta_{max}\}caligraphic_T start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT = { over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ∈ caligraphic_T | italic_δ start_POSTSUBSCRIPT over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT < italic_δ ∧ italic_θ start_POSTSUBSCRIPT over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT < italic_θ start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT } \triangleright Reference set conditions

3:𝒯ref=minΔ(|T~T~|)(T~𝒯ref)(n)\mathcal{T}_{ref}=\displaystyle\min_{\Delta(|\tilde{T}^{\prime}-\tilde{T}|)}(%\tilde{T}^{\prime}\in\mathcal{T}^{\ast}_{ref})_{(n)}caligraphic_T start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT = roman_min start_POSTSUBSCRIPT roman_Δ ( | over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT - over~ start_ARG italic_T end_ARG | ) end_POSTSUBSCRIPT ( over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ∈ caligraphic_T start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT ) start_POSTSUBSCRIPT ( italic_n ) end_POSTSUBSCRIPT \triangleright Reference frame selection

4:forT~𝒯refsuperscript~𝑇subscript𝒯𝑟𝑒𝑓\tilde{T}^{\prime}\in\mathcal{T}_{ref}over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ∈ caligraphic_T start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPTdo

5:T^superscript^𝑇absent\hat{T}^{\prime}\leftarrowover^ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ← SP(T~,f,θminsuperscript~𝑇𝑓subscript𝜃𝑚𝑖𝑛\tilde{T}^{\prime},f,\theta_{min}over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_f , italic_θ start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT) \triangleright Preprocess via spatial channel

6:Δμc=𝒑Pc|T~𝒑T~𝒑||Pc|cΔsubscript𝜇𝑐subscript𝒑subscript𝑃𝑐subscript~𝑇𝒑subscriptsuperscript~𝑇𝒑subscript𝑃𝑐for-all𝑐\Delta\mu_{c}=\frac{\sum_{\boldsymbol{p}\in P_{c}}|\tilde{T}_{\boldsymbol{p}}-%\tilde{T}^{\prime}_{\boldsymbol{p}}|}{|P_{c}|}~{}\forall croman_Δ italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = divide start_ARG ∑ start_POSTSUBSCRIPT bold_italic_p ∈ italic_P start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_POSTSUBSCRIPT | over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT - over~ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT | end_ARG start_ARG | italic_P start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT | end_ARG ∀ italic_c

7:T~𝒑T~𝒑+Δμc𝒑Pccsuperscriptsubscript~𝑇𝒑superscriptsubscript~𝑇𝒑Δsubscript𝜇𝑐for-all𝒑subscript𝑃𝑐for-all𝑐\tilde{T}_{\boldsymbol{p}}^{\prime}\leftarrow\tilde{T}_{\boldsymbol{p}}^{%\prime}+\Delta\mu_{c}~{}\forall\boldsymbol{p}\in P_{c}~{}\forall cover~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ← over~ start_ARG italic_T end_ARG start_POSTSUBSCRIPT bold_italic_p end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT + roman_Δ italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ∀ bold_italic_p ∈ italic_P start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ∀ italic_c \triangleright Linear adjustments

8:endfor

9:return T^tempaverage(𝒯ref)superscript^𝑇𝑡𝑒𝑚𝑝averagesubscript𝒯𝑟𝑒𝑓\hat{T}^{temp}\leftarrow\text{average}(\mathcal{T}_{ref})over^ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT italic_t italic_e italic_m italic_p end_POSTSUPERSCRIPT ← average ( caligraphic_T start_POSTSUBSCRIPT italic_r italic_e italic_f end_POSTSUBSCRIPT ) \triangleright Interpolate

We set n=3𝑛3n=3italic_n = 3, δ=2𝛿2\delta=2italic_δ = 2, and θmax=0.1subscript𝜃max0.1\theta_{\text{max}}=0.1italic_θ start_POSTSUBSCRIPT max end_POSTSUBSCRIPT = 0.1. SP(\cdot) denotes the function for spatial channel defined in Algorithm 1, and f𝑓fitalic_f and θminsubscript𝜃min\theta_{\text{min}}italic_θ start_POSTSUBSCRIPT min end_POSTSUBSCRIPT at line 4 follow their default values from Sec.4.1.The subscript (n)subscript𝑛\cdot_{(n)}⋅ start_POSTSUBSCRIPT ( italic_n ) end_POSTSUBSCRIPT used in line 2 signifies the first n𝑛nitalic_n elements of the sequence, arranged according to the minimization criterion.Pcsubscript𝑃𝑐P_{c}italic_P start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is the set of all pixel locations where the corresponding land cover class is c𝑐citalic_c, and |Pc|subscript𝑃𝑐|P_{c}|| italic_P start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT | is the cardinality of set Pcsubscript𝑃𝑐P_{c}italic_P start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT.

4.3 Estimate LST via Weighted Average

Following the previous two subsections, we acquire initial predictions T^spsuperscript^𝑇𝑠𝑝\hat{T}^{sp}over^ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT italic_s italic_p end_POSTSUPERSCRIPT from the spatial channel and T^tempsuperscript^𝑇𝑡𝑒𝑚𝑝\hat{T}^{temp}over^ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT italic_t italic_e italic_m italic_p end_POSTSUPERSCRIPT from the temporal channel. The final interpolated LST, denoted as T^^𝑇\hat{T}over^ start_ARG italic_T end_ARG, is calculated as the weighted average of the two initial predictions:

T^=wT^sp+(1w)T^temp^𝑇𝑤superscript^𝑇𝑠𝑝1𝑤superscript^𝑇𝑡𝑒𝑚𝑝\hat{T}=w\cdot\hat{T}^{sp}+(1-w)\cdot\hat{T}^{temp}over^ start_ARG italic_T end_ARG = italic_w ⋅ over^ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT italic_s italic_p end_POSTSUPERSCRIPT + ( 1 - italic_w ) ⋅ over^ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT italic_t italic_e italic_m italic_p end_POSTSUPERSCRIPT(8)

where we set the weight w=1θ𝑤1𝜃w=1-\thetaitalic_w = 1 - italic_θ.Note that for a minimally occluded image (a small θ𝜃\thetaitalic_θ), more weight is assigned to the spatial channel.Contrarily, for a severely occluded image (a large θ𝜃\thetaitalic_θ), fewer neighboring pixels are available for the spatial channel, and more weight is assigned to the temporal channel accordingly.

5 Results

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

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).333For 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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG every 16 days444Subject to data availability and requires θ<0.99𝜃0.99\theta<0.99italic_θ < 0.99. for each observation region.Our public dataset is available on the NHERI DesignSafeCyberinfrastructure555https://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).

ParametersMAE (K) \downarrowRMSE (K) \downarrow
Citys𝑠sitalic_sn𝑛nitalic_nM1M2M3M4M5M1M2M3M4M5
Houston250101.882.191.822.494.262.432.852.333.105.05
Houston75032.002.391.872.573.972.533.052.383.254.71
Jacksonville7521.471.571.532.143.891.882.071.882.614.51
Phoenix50011.962.081.362.182.632.622.781.812.883.33
New York10021.251.371.142.113.821.651.841.422.634.28

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×s𝑠𝑠s\times sitalic_s × italic_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𝑛nitalic_n occlusion regions of size s×s𝑠𝑠s\times sitalic_s × italic_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.62Ktimes2.62K2.62\text{\,}\mathrm{K}start_ARG 2.62 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG 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 study666Not 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 θ𝜃\thetaitalic_θ), 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.

ISLAND: Interpolating Land Surface Temperature using land cover (5)
ISLAND: Interpolating Land Surface Temperature using land cover (6)

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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG 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.95limit-from2.952.95-2.95 -2.99Ktimes2.99kelvin2.99\text{\,}\mathrm{K}start_ARG 2.99 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG across four sites.777We 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.55Ktimes3.55kelvin3.55\text{\,}\mathrm{K}start_ARG 3.55 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG and 4.78Ktimes4.78kelvin4.78\text{\,}\mathrm{K}start_ARG 4.78 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG 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.90limit-from3.903.90-3.90 -4.75Ktimes4.75kelvin4.75\text{\,}\mathrm{K}start_ARG 4.75 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG, across four sites. Compared to clear-sky RMSE in Fig.5, ISLAND introduces an additional 0.95limit-from0.950.95-0.95 -1.76Ktimes1.76kelvin1.76\text{\,}\mathrm{K}start_ARG 1.76 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG RMSE error across four sites, with the average additional RMSE being 1.33Ktimes1.33kelvin1.33\text{\,}\mathrm{K}start_ARG 1.33 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG.

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.96Ktimes3.96K3.96\text{\,}\mathrm{K}start_ARG 3.96 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG, 5.07Ktimes5.07K5.07\text{\,}\mathrm{K}start_ARG 5.07 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG, 4.88Ktimes4.88K4.88\text{\,}\mathrm{K}start_ARG 4.88 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG, and 3.86Ktimes3.86K3.86\text{\,}\mathrm{K}start_ARG 3.86 end_ARG start_ARG times end_ARG start_ARG roman_K end_ARG, respectively.

Finally, Fig.6 includes reconstruction error for all cloudy-sky conditions except for dates with more than 99% of pixels occluded (θ>0.99𝜃0.99\theta>0.99italic_θ > 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.

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

5.5.1 Deriving surface temperature trends

The robustness demonstrated in Sec.5.25.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 θ<.99𝜃.99\theta<.99italic_θ < .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 (θ<.99𝜃.99\theta<.99italic_θ < .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 15superscript1515^{\circ}15 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPTC 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.

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

5.5.2 Urban heat island effects

Another key application of our model is the ability to study urban heat islandeffects (UHIE)(Moller etal.,, 2022) at high spatiotemporal resolution.UHIE refers to the phenomenon of urban areas being significantly warmer than their surrounding rural areas.UHIE is primarily driven by the differences in thermal absorption between different materials.For example, grass- or water-covered surfaces tend to have lower temperatures than concrete and asphalt.By providing high spatial and temporal resolution LST outputs, ISLAND offers a novel data product for identifying, studying, and monitoring UHIE in major U.S. metropolitan areas.Fig.8 shows maps of three of the largest U.S. metropolitan areas, Los Angeles, Chicago, and Houston, where the pixel values indicate the number of days surpassing a region-specific temperature threshold.The thresholds are selected based on the definition of Extreme Dangerfrom the National Weather Service (NWS) heat index(Rothfusz and Headquarters,, 1990).The NWS heat index is a function of both temperature and relative humidity.The Comparative Climatic Data (CCD-2018)(NOAA,, 2020) providesthe morning annual average relative humidity (R𝑅Ritalic_R) for each city.We select the LST threshold for each region as the minimum ambient dry bulb temperature that meets the NWS Extreme Danger criteria for the city’s corresponding R𝑅Ritalic_R.The range of observations is 4.5years, at a sampling rate of once per 16 days.Higher values in the frequency maps indicate a more frequent occurrence of UHIE.These frequency maps are available at 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG resolution and can be easily computed using our public dataset.From an urban planning perspective, these UHIE frequency maps offer a powerful tool for enhancing our understanding of how land cover choices influence micro-climates, heat extremes, and the associated health risks.By providing insights into the spatial distribution and frequency of UHIE, these maps can inform decision-making processes regarding urban development and land cover management, aiming to mitigate the adverse effects of heat on public health and well-being.

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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG ×\times× 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG)exceeds a land surface temperature threshold of 35superscript3535^{\circ}35 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPTC(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𝑟0.58r=0.58italic_r = 0.58). A similar trend can be seen in cities such as San Antonio (r=0.45𝑟0.45r=0.45italic_r = 0.45) and San Francisco (r=0.46𝑟0.46r=0.46italic_r = 0.46), San Jose (r=0.44𝑟0.44r=0.44italic_r = 0.44), and New York (r=0.34𝑟0.34r=0.34italic_r = 0.34). In contrast, cities such as Houston (r=0.01𝑟0.01r=-0.01italic_r = - 0.01), Dallas (r=0.10𝑟0.10r=-0.10italic_r = - 0.10), Jacksonville (r=0.17𝑟0.17r=0.17italic_r = 0.17), and Fort Worth (r=0.08𝑟0.08r=0.08italic_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.

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

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 Discussion & Conclusions

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 at1kmtimes1km1\text{\,}\mathrm{k}\mathrm{m}start_ARG 1 end_ARG start_ARG times end_ARG start_ARG roman_km end_ARG).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

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

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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG, 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 1kmtimes1kilometer1\text{\,}\mathrm{km}start_ARG 1 end_ARG start_ARG times end_ARG start_ARG roman_km end_ARG.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 7792mtimes7792meter7792\text{\,}\mathrm{m}start_ARG 7792 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG.

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 3kmtimes3kilometer3\text{\,}\mathrm{km}start_ARG 3 end_ARG start_ARG times end_ARG start_ARG roman_km end_ARG andERA5(MuñozSabater,, 2019) at around 11kmtimes11kilometer11\text{\,}\mathrm{km}start_ARG 11 end_ARG start_ARG times end_ARG start_ARG roman_km end_ARG.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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG 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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG.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 30mtimes30meter30\text{\,}\mathrm{m}start_ARG 30 end_ARG start_ARG times end_ARG start_ARG roman_m end_ARG 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.

Acknowledgement

The authors gratefully acknowledge the support of this research by the National Science Foundation (NSF) award numbers 1652633 and 2107313. The contributions of Pranavesh Panakkal and Jamie E. Padgett were partially supported by NSF award number 2227467. Any opinions, findings, conclusions, or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the sponsors.

Appendix A Additional Details on NLCD Land Cover

Fig.11 shows the visualization of the NLCD maps referenced in this paper. Fig.11(a) shows urban NLCD maps for cities in our visual results (Fig.5.1), simulation evaluation (Sections5.25.3), and illustrated applications (Fig.5.5), while Fig.11(b) shows the NLCD maps around four SURFRAD sites for our in situ evaluation (Fig.5.4). Refer to Fig.11(c) for the legend of NLCD classes.

The heterogeneity of land cover types in urban settings is clearly reflected in Fig.11(a). As discussed in Fig.2 and 3, LST is closely related to land cover; as such, complex terrains in cities make LST reconstruction challenging. The difference in land cover distributions between urban regions is also evident in Fig.11(a). New York City is characterized by high-density developments with very sparse natural landscapes. Houston, TX, has extensive urban and suburban developments with mixed land use, accompanied by agricultural land, forest, and water bodies on the outskirts. Jacksonville, FL, has large expanses of coastal wetlands. Phoenix, AZ, is an urban region surrounded by deserts and mountains with sparse vegetation. We chose the four regions to conduct our simulation evaluation (Table1) to represent ISLAND’s performance under a variety of settings.

Our in situ validation targets are SURFRAD stations located in rural regions. Fig.11(b) shows the relative hom*ogeneous land cover types surrounding the SURFRAD stations, which is in stark contrast to the urban regions in Fig.11(a).

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

References

  • (1)Augustine, J.A., DeLuisi, J.J., and Long, C.N. (2000a).Surfrad - a national surface radiation budget network for atmospheric research.Bulletin of the American Meteorological Society, 81(10):2341–2357.
  • (2)Augustine, J.A., DeLuisi, J.J., and Long, C.N. (2000b).Surfrad–a national surface radiation budget network for atmospheric research.Bulletin of the American Meteorological Society, 81(10):2341–2358.
  • Baiocchi etal., (2017)Baiocchi, V., Zottele, F., and Dominici, D. (2017).Remote sensing of urban microclimate change in l’aquila city (italy) after post-earthquake depopulation in an open source GIS environment.Sensors, 17.
  • Buchhorn etal., (2020)Buchhorn, M., Lesiv, M., Tsendbazar, N.-E., Herold, M., Bertels, L., and Smets, B. (2020).Copernicus global land cover layers—collection 2.Remote Sensing, 12(6):1044.
  • Centers for Disease Control etal., (2020)Centers for Disease Control, Prevention and ATSDR, Analysis and Services Program, and Agency for Toxic Substances and Disease Registry/ Geospatial Research (2020).CDC/ATSDR Social Vulnerability Index 2020 Database US.
  • Chaudhuri and Mishra, (2016)Chaudhuri, G. and Mishra, N.B. (2016).Spatio-temporal dynamics of land cover and land surface temperature in ganges-brahmaputra delta: A comparative analysis between india and bangladesh.Applied Geography, 68:68–83.
  • Chun and Guldmann, (2018)Chun, B. and Guldmann, J.-M. (2018).Impact of greening on the urban heat island: Seasonal variations and mitigation strategies.Computers, Environment and Urban Systems, 71:165–176.
  • Dee etal., (2022)Dee, S.G., Nabizadeh, E., Nittrouer, C.L., Baldwin, J.W., Li, C., Gaviria, L., Guo, S., Lu, K., Saunders-Shultz, B.M., Gurwitz, E., etal. (2022).Increasing health risks during outdoor sports due to climate change in texas: Projections versus attitudes.GeoHealth, 6(8):e2022GH000595.
  • Demarquet etal., (2023)Demarquet, Q., Rapinel, S., Dufour, S., and Hubert-Moy, L. (2023).Long-term wetland monitoring using the landsat archive: A review.Remote Sensing, 15(3):820.
  • Dewitz, (2021)Dewitz, J. (2021).National land cover database (NLCD) 2019 products.
  • Donlon etal., (2012)Donlon, C., Berruti, B., Buongiorno, A., Ferreira, M.-H., Féménias, P., Frerick, J., Goryl, P., Klein, U., Laur, H., Mavrocordatos, C., etal. (2012).The global monitoring for environment and security (gmes) sentinel-3 mission.Remote sensing of Environment, 120:37–57.
  • Dugord etal., (2014)Dugord, P.-A., Lauf, S., Schuster, C., and Kleinschmit, B. (2014).Land use patterns, temperature distribution, and potential heat stress risk–the case study berlin, germany.Computers, Environment and Urban Systems, 48:86–98.
  • Eisavi etal., (2016)Eisavi, V., Yazdi, A.M., and Niknezhad, S.A. (2016).Spatial and temporal modeling of wetland surface temperature using landsat-8 imageries in sulduz, iran.Journal of the Faculty of Forestry Istanbul University, 66(1):46–58.
  • Fayad etal., (2021)Fayad, F.A., Maref, W., and Awad, M.M. (2021).Review of white roofing materials and emerging economies with focus on energy performance cost-benefit, maintenance, and consumer indifference.Sustainability, 13(17):9967.
  • Gomez-Martinez etal., (2021)Gomez-Martinez, F., deBeurs, K.M., Koch, J., and Widener, J. (2021).Multi-temporal land surface temperature and vegetation greenness in urban green spaces of puebla, mexico.Land, 10(2):155.
  • Gorelick etal., (2017)Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., and Moore, R. (2017).Google earth engine: Planetary-scale geospatial analysis for everyone.Remote Sensing of Environment.
  • Hakim etal., (2016)Hakim, G.J., Emile-Geay, J., Steig, E.J., Noone, D., Anderson, D.M., Tardif, R., Steiger, N., and Perkins, W.A. (2016).The last millennium climate reanalysis project: Framework and first results.Journal of Geophysical Research: Atmospheres, 121(12):6745–6764.
  • Han etal., (2024)Han, J., Fang, S., Mi, Q., Wang, X., Yu, Y., Zhuo, W., and Peng, X. (2024).A time-continuous land surface temperature (lst) data fusion approach based on deep learning with microwave remote sensing and high-density ground truth observations.Science of The Total Environment, 914:169992.
  • Horton etal., (2016)Horton, R.M., Mankin, J.S., Lesk, C., Coffel, E., and Raymond, C. (2016).A review of recent advances in research on extreme heat events.Current Climate Change Reports, 2:242–259.
  • Huang and Wang, (2019)Huang, X. and Wang, Y. (2019).Investigating the effects of 3d urban morphology on the surface urban heat island effect in urban functional zones by using high-resolution remote sensing data: A case study of wuhan, central china.ISPRS Journal of Photogrammetry and Remote Sensing, 152:119–131.
  • Hulley etal., (2015)Hulley, G.C., Hook, S.J., Abbott, E., Malakar, N., Islam, T., and Abrams, M. (2015).The ASTER global emissivity dataset (ASTER GED): Mapping earth’s emissivity at 100 meter spatial scale.Geophysical Research Letters, 42(19):7966–7976.
  • Imran etal., (2021)Imran, H.M., Hossain, A., Islam, A. K. M.S., Rahman, A., Bhuiyan, M. A.E., Paul, S., and Alam, A. (2021).Impact of land cover changes on land surface temperature and human thermal comfort in dhaka city of bangladesh.Earth Systems and Environment, 5(3):667–693.
  • James etal., (2022)James, E.P., Alexander, C.R., Dowell, D.C., Weygandt, S.S., Benjamin, S.G., Manikin, G.S., Brown, J.M., Olson, J.B., Hu, M., Smirnova, T.G., etal. (2022).The high-resolution rapid refresh (hrrr): An hourly updating convection-allowing forecast model. part ii: Forecast performance.Weather and Forecasting, 37(8):1397–1417.
  • Jin etal., (2019)Jin, Z., Zhang, Y., DelGenio, A., Schmidt, G., and Kelley, M. (2019).Cloud scattering impact on thermal radiative transfer and global longwave radiation.Journal of quantitative spectroscopy & radiative transfer, 239:106669.
  • Knapp etal., (2011)Knapp, K.R., Ansari, S., Bain, C.L., Bourassa, M.A., Dickinson, M.J., Funk, C., Helms, C.N., Hennon, C.C., Holmes, C.D., Huffman, G.J., Kossin, J.P., Lee, H.-T., Loew, A., and Magnusdottir, G. (2011).Globally gridded satellite observations for climate studies.Bulletin of the American Meteorological Society, 92(7):893–907.
  • Knapp, (2014)Knapp, KennethR., N. C.P. (2014).Noaa climate data record (cdr) of intersatellite calibrated gridded satellite data from isccp b1 (gridsat-b1) 11 micron brightness temperature, version 2.
  • Li etal., (2013)Li, Z.-L., Tang, B.-H., Wu, H., Ren, H., Yan, G., Wan, Z., Trigo, I.F., and Sobrino, J.A. (2013).Satellite-derived land surface temperature: Current status and perspectives.Remote Sensing of Environment, 131:14–37.
  • Malakar etal., (2018)Malakar, N.K., Hulley, G.C., Hook, S.J., Laraby, K., Cook, M., and Schott, J.R. (2018).An operational land surface temperature product for landsat thermal data: Methodology and validation.IEEE Transactions on Geoscience and Remote Sensing, 56(10):5717–5735.
  • McShane etal., (2022)McShane, C., Uhl, J.H., and Leyk, S. (2022).Gridded land use data for the conterminous united states 1940–2015.Scientific Data, 9(1):493.
  • Moller etal., (2022)Moller, V., van Diemen, R., Matthews, J., Méndez, C., sem*nov, S., Fuglestvedt, J., and Reisinger, A. (2022).Annex ii: Glossary.In Pörtner, H.-O., Roberts, D., Tignor, M., Poloczanska, E., Mintenbeck, K., Alegría, A., Craig, M., Langsdorf, S., Löschke, S., Möller, V., Okem, A., and Rama, B., editors, Climate Change 2022: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, pages 2897–2930. Cambridge University Press, Cambridge, UK and New York, NY, USA.
  • MuñozSabater, (2019)MuñozSabater, J. (2019).Era5-land monthly averaged data from 2001 to present.
  • NOAA, (2020)NOAA (2020).Comparative climatic data.
  • Ogawa etal., (2008)Ogawa, K., Schmugge, T., and Rokugawa, S. (2008).Estimating broadband emissivity of arid regions and its seasonal variations using thermal infrared remote sensing.IEEE Transactions on Geoscience and Remote Sensing, 46(2):334–343.
  • Orimoloye etal., (2018)Orimoloye, I.R., Mazinyo, S.P., Nel, W., and Kalumba, A.M. (2018).Spatiotemporal monitoring of land surface temperature and estimated radiation using remote sensing: human health implications for east london, south africa.Environmental Earth Sciences, 77(3):77.
  • Osborne and Alvares-Sanches, (2019)Osborne, P.E. and Alvares-Sanches, T. (2019).Quantifying how landscape composition and configuration affect urban land surface temperatures using machine learning and neutral landscapes.Computers, Environment and Urban Systems, 76:80–90.
  • Paris etal., (2009)Paris, S., Kornprobst, P., Tumblin, J., Durand, F., etal. (2009).Bilateral filtering: Theory and applications.Foundations and Trends® in Computer Graphics and Vision, 4(1):1–73.
  • Qiu etal., (2020)Qiu, S., Zhu, Z., and Woodco*ck, C.E. (2020).Cirrus clouds that adversely affect landsat 8 images: What are they and how to detect them?Remote Sensing of Environment, 246:111884.
  • Rathje etal., (2017)Rathje, E.M., Dawson, C., Padgett, J.E., Pinelli, J.-P., Stanzione, D., Adair, A., Arduino, P., Brandenberg, S.J., co*ckerill, T., Dey, C., etal. (2017).Designsafe: New cyberinfrastructure for natural hazards engineering.Natural Hazards Review, 18(3):06017001.
  • Rothfusz and Headquarters, (1990)Rothfusz, L.P. and Headquarters, N. S.R. (1990).The heat index equation (or, more than you ever wanted to know about heat index).Fort Worth, Texas: National Oceanic and Atmospheric Administration, National Weather Service, Office of Meteorology, 9023:640.
  • Roy etal., (2008)Roy, D.P., Ju, J., Lewis, P., Schaaf, C., Gao, F., Hansen, M., and Lindquist, E. (2008).Multi-temporal modis–landsat data fusion for relative radiometric normalization, gap filling, and prediction of landsat data.Remote Sensing of Environment, 112(6):3112–3130.
  • Seneviratne etal., (2021)Seneviratne, S., Zhang, X., Adnan, M., Badi, W., Dereczynski, C., DiLuca, A., Ghosh, S., Iskandar, I., Kossin, J., Lewis, S., Otto, F., Pinto, I., Satoh, M., Vicente-Serrano, S., Wehner, M., and Zhou, B. (2021).Weather and Climate Extreme Events in a Changing Climate, page 1513–1766.Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA.
  • Sobrino etal., (2013)Sobrino, J.A., Oltra-Carrió, R., Sòria, G., Jiménez-Muñoz, J.C., Franch, B., Hidalgo, V., Mattar, C., Julien, Y., Cuenca, J., Romaguera, M., Gómez, J.A., DeMiguel, E., Bianchi, R., and Paganini, M. (2013).Evaluation of the surface urban heat island effect in the city of madrid by thermal remote sensing.International Journal of Remote Sensing, 34(9):3177–3192.Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/01431161.2012.716548.
  • Stillinger etal., (2019)Stillinger, T., Roberts, D.A., Collar, N.M., and Dozier, J. (2019).Cloud masking for landsat 8 and modis terra over snow-covered terrain: Error analysis and spectral similarity between snow and cloud.Water Resources Research, 55(7):6169–6184.
  • Streutker, (2003)Streutker, D.R. (2003).Satellite-measured growth of the urban heat island of houston, texas.Remote Sensing of Environment, 85(3):282–289.
  • The White House, (2022)The White House (2022).Justice40 Initiative | Environmental Justice | The White House — whitehouse.gov.https://www.whitehouse.gov/environmentaljustice/justice40/.[Accessed 22-May-2023].
  • Tobler, (1970)Tobler, W.R. (1970).A computer movie simulating urban growth in the detroit region.Economic Geography, 46:234–240.
  • United States Census Bureau, (2023)United States Census Bureau (2023).City and town population totals: 2020-2021.
  • (48)USGS (2023a).Landsat-8 image courtesy of the u.s. geological survey.Accessed: April 3, 2023.
  • (49)USGS (2023b).Landsat-9 image courtesy of the u.s. geological survey.Accessed: April 3, 2023.
  • Vo etal., (2023)Vo, T.T., Hu, L., Xue, L., Li, Q., and Chen, S. (2023).Urban effects on local cloud patterns.Proceedings of the National Academy of Sciences, 120(21):e2216765120.
  • Wan, (2013)Wan, Z. (2013).Modis land surface temperature products.
  • Wang and Liang, (2009)Wang, K. and Liang, S. (2009).Evaluation of aster and modis land surface temperature and emissivity products using long-term surface longwave radiation observations at surfrad sites.Remote Sensing of Environment, 113(7):1556–1565.
  • Wang etal., (2019)Wang, T., Shi, J., Ma, Y., Husi, L., Comyn-Platt, E., Ji, D., Zhao, T., and Xiong, C. (2019).Recovering land surface temperature under cloudy skies considering the solar-cloud-satellite geometry: Application to modis and landsat-8 data.Journal of Geophysical Research: Atmospheres, 124(6):3401–3416.
  • Weng and Fu, (2014)Weng, Q. and Fu, P. (2014).Modeling annual parameters of clear-sky land surface temperature variations and evaluating the impact of cloud cover using time series of landsat tir data.Remote Sensing of Environment, 140:267–278.
  • Wu etal., (2021)Wu, P., Yin, Z., Zeng, C., Duan, S.-B., Göttsche, F.-M., Ma, X., Li, X., Yang, H., and Shen, H. (2021).Spatially continuous and high-resolution land surface temperature product generation: A review of reconstruction and spatiotemporal fusion techniques.IEEE Geoscience and Remote Sensing Magazine, 9(3):112–137.
  • Wu, (2020)Wu, Q. (2020).geemap: A python package for interactive mapping with google earth engine.Journal of Open Source Software, 5(51):2305.
  • Yang and Huang, (2021)Yang, J. and Huang, X. (2021).The 30m annual land cover dataset and its dynamics in china from 1990 to 2019.Earth System Science Data, 13(8):3907–3925.
  • Yang etal., (2018)Yang, L., Jin, S., Danielson, P., Homer, C., Gass, L., Bender, S.M., Case, A., Costello, C., Dewitz, J., Fry, J., Funk, M., Granneman, B., Liknes, G.C., Rigge, M., and Xian, G. (2018).A new generation of the united states national land cover database: Requirements, research priorities, design, and implementation strategies.ISPRS Journal of Photogrammetry and Remote Sensing, 146:108–123.
  • Yu etal., (2019)Yu, W., Tan, J., Ma, M., Li, X., She, X., and Song, Z. (2019).An effective similar-pixel reconstruction of the high-frequency cloud-covered areas of southwest china.Remote Sensing, 11(3).
  • Zeng etal., (2015)Zeng, C., Shen, H., Zhong, M., Zhang, L., and Wu, P. (2015).Reconstructing modis lst based on multitemporal classification and robust regression.IEEE Geoscience and Remote Sensing Letters, 12(3):512–516.
  • Zhao etal., (2020)Zhao, J., Zhao, X., Liang, S., Zhou, T., Du, X., Xu, P., and Wu, D. (2020).Assessing the thermal contributions of urban land cover types.Landscape and Urban Planning, 204:103927.
  • Zhu etal., (2022)Zhu, X., Duan, S.-B., Li, Z.-L., Wu, P., Wu, H., Zhao, W., and Qian, Y. (2022).Reconstruction of land surface temperature under cloudy conditions from landsat 8 data using annual temperature cycle model.Remote Sensing of Environment, 281:113261.
  • Zhu and Woodco*ck, (2012)Zhu, Z. and Woodco*ck, C.E. (2012).Object-based cloud and cloud shadow detection in landsat imagery.Remote sensing of environment, 118:83–94.
ISLAND: Interpolating Land Surface Temperature using land cover (2024)
Top Articles
Latest Posts
Article information

Author: Gregorio Kreiger

Last Updated:

Views: 6520

Rating: 4.7 / 5 (77 voted)

Reviews: 92% of readers found this page helpful

Author information

Name: Gregorio Kreiger

Birthday: 1994-12-18

Address: 89212 Tracey Ramp, Sunside, MT 08453-0951

Phone: +9014805370218

Job: Customer Designer

Hobby: Mountain biking, Orienteering, Hiking, Sewing, Backpacking, Mushroom hunting, Backpacking

Introduction: My name is Gregorio Kreiger, I am a tender, brainy, enthusiastic, combative, agreeable, gentle, gentle person who loves writing and wants to share my knowledge and understanding with you.