See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/325231273

An Operational Land Surface Temperature Product for Landsat Thermal Data:

Methodology and Validation

Article  in  IEEE Transactions on Geoscience and Remote Sensing · May 2018

DOI: 10.1109/TGRS.2018.2824828

CITATIONS READS
9 1,338

6 authors, including:

Nabin K. Malakar Kelly G Laraby

Worcester State University Rochester Institute of Technology
32 PUBLICATIONS   213 CITATIONS    2 PUBLICATIONS   14 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

NASA MODIS Land Surface Temperature and Emissivity Product (MOD21) View project

Special Issue "Sustainability in the Mountains Region" View project

All content following this page was uploaded by Nabin K. Malakar on 31 May 2018.

The user has requested enhancement of the downloaded file.

 An Operational Land Surface Temperature Product for Landsat Thermal Data:
Methodology and Validation

Nabin Malakar, Glynn Hulley, Simon Hook

Jet Propulsion Laboratory, Caltech, Pasadena CA 91109
Kelly Laraby, Monica Cook, John Schott
Rochester Institute of Technology, Rochester, NY 14623

Abstract

Thermal sensors onboard Landsat satellites have been underutilized due to lack of consistent and

accurate methodologies for retrieving the Land Surface Temperature (LST) at global scales over

all land cover types. We present an operational algorithm for generating Landsat LST

consistently for all sensors that will be implemented by the USGS/NASA and made available at

the LPDAAC. The LST algorithm involves three steps: the observed thermal radiance is

atmospherically corrected using a radiative transfer model and reanalysis data; the ASTER

Global Emissivity Dataset v3 (GEDv3) is spectrally adjusted, and then modified to account for

vegetation phenology and snow cover using Landsat VSWIR data; LST is retrieved by inverting

the atmospherically and emissivity corrected Landsat radiances with a look-up-table approach.

Landsat derived emissivities were validated at two pseudo invariant sand dune sites within an

average absolute error of 0.54% when compared to lab measurements. The Landsat LST

retrievals were validated with in-situ observations from four SURFRAD sites, and two inland

water bodies (Salton Sea, Lake Tahoe) in the USA. The LST retrievals for Landsat 5 and Landsat

7 had a mean bias (RMSE) of 0.7 K (2.2 K), and 0.9 K (2.3 K) for the SURFRAD sites, and −0.3

K (0.6 K) and 0.4 K (0.7 K) for the inland water bodies, respectively. The operational algorithm

will provide a consistent LST record from four decades of historical Landsat thermal data

enabling the long-term monitoring of temperature and trends, land cover and land use changes,

and improved utilization in models.

1
1 Introduction
Land surface temperature (LST) is a key environmental climate variable derived from

Thermal Infrared (TIR) data that is used in surface energy balance models in studies involving

the hydrological cycle [1-3], evapotranspiration, drought, and climate research [4-8]. LST is also

useful for the study of heat-related issues, impacts of heat stress on the urban population [9, 10],

and outbreak and propagation of vector-borne diseases [11, 12]. NASA and many other

organizations have identified LST as an important Earth Surface Data Record (ESDR) [13-15].

More recently, with support from the International Land Surface Temperature and Emissivity

Working Group (ILSTE-WG, http://ilste-wg.org/), LST was accepted and defined as an

Environmental Climate Variable (ECV) by the Global Climate Observing System (GCOS).

The Landsat suite of satellites has been one of the major contributing factors in the

development of global-scale earth systems science research such as the international geosphere-

biosphere program (IGBP) [16]. In particular, data from the Landsat thermal bands beginning

with Landsat 4 in 1982, presents a unique opportunity to study long-term trends of Earth surface

temperatures at ~100m spatial resolution that could be used for monitoring climate change at

regional and local (urban) scales. Table 1 presents the Landsat suite of TIR measurement

characteristics. Landsat 6 failed to reach the orbit, and has been omitted in the discussions. The

Landsat 4 and 5 satellites carried both the Multispectral Scanner (MSS) and the Thematic

Mapper (TM) sensors. After the decommissioning of Landsat 5 in 2013, two satellites are

currently in operation: Landsat 7 and Landsat 8. The Landsat 7 satellite was launched in 1999

and carries the Enhanced Thematic Mapper Plus (ETM+) sensor which has a spatial resolution of

60m in the TIR. The ETM+ Scan Line Corrector (SLC) failed in May 2003, but still continues to

gather data. The SLC failure has been reported to not influence the calibration of the thermal

2
band [17]. The Landsat 8 satellite was launched in 2013 and carries the Operational Land Imager

(OLI) and the Thermal Infrared Sensor (TIRS). The two TIR bands of Landsat 8 allow a split

window approach; however, an offset problem with Landsat 8 calibrated TIR bands was

identified by USGS and NASA in August 2013. The observed offset is introduced by stray light

entering the TIRS sensors field of view [18]. Due to a large calibration uncertainty associated

with Landsat 8 TIR band 11, USGS recommended that it should not be used for split window

based retrievals of LST [19].

The moderate-level spatial resolution offered by the Landsat TIR bands can sufficiently

capture land cover changes at local scales, e.g. due to anthropogenic and natural factors [20].

However, the thermal infrared data from Landsat have been largely under-utilized due to the lack

of an accurate, standard Landsat LST product. The difficulty in producing such a product lies

mostly in the fact that most Landsat thermal sensors have only one single thermal band (Landsat

4-7). The radiance observed in the TIR is a function of both the temperature and emissivity of the

surface, and since the Landsat 4-7 sensors have only a single thermal band, the problem of

separating the contributions from temperature and emissivity from the observed radiance

becomes an ill-posed problem (i.e. two variables to solve for – temperature and emissivity –

from one piece of information – the observed thermal radiance). Past studies [21-23] derived

LST by assigning surface emissivity based on land cover classification schemes and vegetation

indices [24], but these approaches are most suited for local studies and can be problematic when

applied globally, in particular over arid and semi-arid regions where the emissivities vary widely,

both spectrally and spatially. For example, an emissivity error of 0.015 (1.5%) will result in an

LST error of ~1 K for a temperature of 300 K at 11 µm [25].

3
 In 2003, Barsi et al. [26] developed a web-based tool to derive atmospheric parameters

using input National Centers for Environmental Prediction (NCEP) atmospheric profile and

commercially available MODTRAN software. The interface provided the site specific

atmospheric transmission, upwelling and downwelling radiances for the targets with known

emissivity and clear sky conditions. These parameters could be applied to atmospherically

correct single thermal band of Landsat 5, Landsat 7 and Landsat 8 (band 10) to obtain the LST.

Validation of the LST with in-situ data for selected dates showed that the retrieval was accurate

within a fraction of degrees for the selected scenes under clear sky conditions [17]. However,

the web-based method had some limitations; atmospheric parameters were generated for a single

NCEP 1°×1° grid point closest to the desired scene, the NCEP record was not available for the

entire Landsat lifetime, and the user had to apply their own emissivity correction [17] to estimate

the LST.

In 2014, a new source of emissivity information became available with the release of the

Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global

Emissivity Dataset (GED) version 3 (v3) [27]. The ASTER GEDv3 is a global emissivity dataset

produced from using millions of cloud free ASTER emissivity data from 2000-2008. These data

were gridded, averaged, and mosaicked to produce a global dataset for the five ASTER TIR

wavelengths (8.3, 8.6, 9.1, 10.6, and 11.3 µm) and for two different resolutions - 3 arc sec (~100

m) and 30 arc sec (~1 km). ASTER GEDv3 has been extensively validated in the past over

mostly arid and semi-arid regions [27, 28] with an average absolute band error of ~1% [27]. The

availability of ASTER GEDv3 at appropriate spatial resolutions for Landsat was a major step

forward in the generation of a standard and consistent Landsat LST product over the entire globe.

The ASTER GEDv3 is first spectrally adjusted to the Landsat TIR bands, and then adjusted for

4
vegetation phenology and seasonal changes in snow cover using information from the Landsat

visible and shortwave infrared data.

The paper is organized as follows: In the section 2, we will describe the basic radiative

transfer theory in the thermal infrared domain, followed by the atmospheric correction method

(section 3), the emissivity adjustments of ASTER GEDv3 to Landsat TIR bands (section 4), and

the LST retrieval method (section 5). The study concludes with a validation study over two

different cover types; one over two water bodies: Lake Tahoe and Salton Sea located in

California/Nevada, USA (section 6), and the other over four SURFRAD land sites in the USA. A

more comprehensive validation of the LST product over a larger variety of land cover types will

be performed in a future study.

2 Theoretical Background
The Earth-emitted radiance is a function of temperature and emissivity, and gets attenuated

by the atmosphere on its path to the satellite’s sensor. The clear-sky top of atmosphere radiance

measured by a spaceborne sensor (Lsat,λ), in the narrow band of the sensor channel defined

between [λ1, λ2] centered on wavelength λ, includes contributions from surface emission,

atmospheric path radiance (L↑λ) and atmospheric downwelling irradiance (L↓λ) reflected by the

Earth’s surface and attenuated by the atmosphere as follows:

𝐿𝐿𝑠𝑠𝑠𝑠𝑠𝑠,𝜆𝜆 (𝜃𝜃) = �𝜀𝜀𝜆𝜆 𝐵𝐵𝜆𝜆 (𝑇𝑇𝑠𝑠 ) + (1 − 𝜀𝜀𝜆𝜆 ) 𝐿𝐿↓𝜆𝜆 � 𝜏𝜏𝜆𝜆 (𝜃𝜃) + 𝐿𝐿↑𝜆𝜆 (𝜃𝜃) (1)

where Bλ(T) is the Planck function describing the radiance of a blackbody at temperature 𝑇𝑇𝑠𝑠 , 𝜖𝜖𝜆𝜆

is effective emissivity and τλ is atmospheric transmissivity, and 𝜃𝜃 is the viewing angle of the

sensor. In Eq. (1), the terms in square brackets represent the total radiance leaving the surface,

which is the sum of direct surface emission and reflected sky radiance. In the longwave TIR, the

reflected solar radiation is much smaller, and negligible. We also assumed that the surface is
5
Lambertian (angular variation in emissivity is small) and used Kirchhoff’s law to express the

hemispherical directional reflectance as directional emissivity (r = 1-.𝜀𝜀𝜆𝜆 ) [29]. The radiance

measurements are typically performed over a discrete sensor band and represent the convolution

of at-sensor spectral radiance with the sensor's spectral response function. Since here only a single

TIR band is used, the subscript λ will be dropped henceforth.

For each Landsat instrument, the sensor spectral response, orbital location, and line of

sight vector for each pixel are known, so the problem reduces to characterizing and

compensating the atmospheric effects. The atmospheric parameters, L↑, L↓ and τ are estimated

with a radiative transfer model using input atmospheric fields of air temperature, relative

humidity, and geopotential height at the time and location of the measurement. However, even

after atmospheric correction, Eq. (1) is underdetermined in that the number of variables to solve

for (one temperature and n emissivities) is always greater than the number of observations

available (n) [30, 31]. As such, multiple retrieval algorithms have been developed for separating

the temperature and emissivity components from the observed radiance using various constraints

for different sensors and configurations [29, 32-36]. In the case of multispectral TIR sensors,

two primary algorithms have been used for estimating LST and emissivity — day-night

method[29] and Temperature and Emissivity Separation (TES) algorithm [34]. The split-

window algorithms assume a fixed emissivity a priori based on land cover classification data and

the differential absorption between two TIR bands for atmospheric correction [37].

Alternatively, the TES algorithm uses full radiative transfer calculations for the atmospheric

correction and exploits an empirical relationship between minimum emissivity and spectral

variance in observed radiance within 3 or more bands [34]. However, since the Landsat satellites

have only one thermal band, the emissivity of the thermal bands has to be derived externally

6
before inverting atmospherically corrected observed radiance to derive the surface temperature

[21, 35, 40-42]. This method is described in the following three sections.

3 Atmospheric Correction
One of the most common models used to simulate atmospheric propagation of

electromagnetic radiation is the MODerate resolution atmospheric TRANsmission (MODTRAN)

radiative transfer model developed by the U.S. Air Force [43-45]. The model requires input

atmospheric profiles of temperature, relative humidity, and geopotential height in addition to

elevation information. MODTRAN has already been used for Landsat calibration studies [46].

For Landsat data over the North American region, we used atmospheric profiles obtained

from the National Weather Service's NCEP reanalysis datasets Global Data Assimilation System

(GDAS) known as North American Regional Reanalysis (NARR) dataset. The 3-hourly 32-km

horizontal resolution archive data from NCEP NARR has 29 pressure levels over the North

American domain and spans the time of 1979 to the present [47]. NARR was primarily

developed to improve on NCEP reanalysis by assimilating precipitation accuracy [48, 49].

Temperature, geopotential height, relative humidity, and ozone profiles are retrieved for each

geopotential height level and precipitable water column values are also read from the dataset.

For Landsat data over the rest of the world, the NASA Modern Era Reanalysis for

Research and Applications Version-2 (MERRA-2) can be used. MERRA-2 provides atmospheric

profiles of temperature, water vapor, and geopotential height at 0.5 x 0.625 degree resolution in

latitude and longitude, respectively, globally for every 6 hours at 42 pressure levels. The

MERRA-2 dataset spans the time domain of 1980-present and is focused on the modern satellite

era using a new version of the Goddard Earth Observing Data Assimilation System Version 5

(GEOS-5) produced by NASA GSF Global Modeling and Assimilation Office (GMAO) [50, 51].
7
MERRA-2 has many core features of the GEOS-5 [52] and includes many of NASA’s Earth

Observation System (EOS), radiosondes, dropsondes, aircraft and in-situ rain rates data among

others.

In order to estimate surface radiance, atmospheric effects first need to be removed from

observed radiance before these can be used to separate emissivity and temperature. It is

important to note that in the longwave infrared regime, atmospheric parameters such as

transmission, upwelling radiance, and downwelling radiance are characteristics of the

atmosphere, while surface emissivity and surface temperature are independent variables.

Assuming 𝜀𝜀 = 1, the observed radiance given by Eq. (1) can be expressed as:

𝐿𝐿𝑠𝑠𝑠𝑠𝑠𝑠 = 𝐵𝐵(𝑇𝑇) 𝜏𝜏 + 𝐿𝐿↑ (2)

which is a simple linear relationship between surface emitted radiance and observed at-sensor

radiance with slope equal to transmittance and intercept equal to upwelling sky radiance. Then,

using two temperatures spanning a wide range, say T1 = 273 K and T2 = 310 K, we could

determine the form of Eq. (2) using simple linear regression, i.e. the transmittance and path

radiance can be determined with two MODTRAN simulations. The downwelling radiance can

then be obtained by solving Eq. (1) for 𝜀𝜀 < 1.

𝐿𝐿𝑠𝑠𝑠𝑠𝑠𝑠 − 𝐿𝐿↑
� � − 𝜀𝜀 𝐵𝐵(𝑇𝑇)
𝜏𝜏
𝐿𝐿↓ = � � (3)
(1 − 𝜀𝜀 )

The emissivity is set to 𝜀𝜀 = 0.9 in the third run of MODTRAN to obtain the last of the three

atmospheric parameters. Note that the temperatures and emissivities assumed in the simulations

are only used for estimating the three band-effective atmospheric parameters. MODTRAN 5.2

was used with a spectral resolution of 1 cm-1 to estimate path radiance, transmissivity and

downward sky radiance that are converted to Landsat band-averaged values by convolving them

8
with the respective sensor spectral response functions. Since generating MODTRAN

atmospheric parameters for every pixel on the scene would be computationally intensive, we

computed these parameters at three ground altitudes for surface temperatures of 273 K, 310 K

and equal to air temperature at land-atmosphere interface by setting T=’000’ for the third

MODTRAN run.

Atmospheric profiles of temperature, geopotential height, pressure, relative humidity and

ozone at given pressure levels are used as input to MODTRAN for the atmospheric correction

step. The digital ground elevation data were obtained from GTOPO30 [53] at 30-arc second

resolution. The NARR data are available on a fixed grid and need to first be subset for a given

Landsat scene dimension. Since the NARR dataset is available at 3-hour time-steps, a time-

period bracketing the Landsat overpass time is subset, from which the atmospheric variables are

linearly interpolated to match Landsat overpass time. The radiative transfer parameters are then

computed at each of the model grid points. Since it is computationally intensive to run

MODTRAN for all Landsat points, the atmospheric parameters are interpolated over the entire

Landsat scene using Shepard’s method [54].

4 Emissivity Correction
Since Landsat sensors have only one thermal band (except for Landsat 8), multispectral

retrieval approaches cannot be applied to Landsat thermal data. Therefore, in the past the primary

limiting factor for generating a consistent Landsat product was the lack of physically derived

emissivity product available at Landsat spatial resolution. Previous approaches have used fixed

emissivities based on land cover classification methods or vegetation threshold methods, and

adjusted these to vegetation dynamics [24, 55, 56]. However, at global scales these approaches

9
are limited since they usually assign a single emissivity to all bare soils and suffer from errors in

land classification products used to assign the emissivity values.

Providing appropriate and accurate background emissivities for LST retrieval has been

greatly improved with the completion of the ASTER GEDv3, which contains the multi-annual

mean gridded global emissivities for all ASTER clear-sky scenes acquired between 2000 and

2008. The product uses a TES algorithm and Water Vapor Scaling (WVS) atmospheric

correction method at ~100 m spatial resolution for all five ASTER TIR bands in the 8-12 µm

region and is available at the Land Processes Distributed Active Archive Center (LP-DAAC)

(https://lpdaac.usgs.gov/products/community_products_table). The accuracy of the ASTER

GEDv3 emissivity is on average ~0.01 for all bands over arid and semi-arid regions [57].

For use with Landsat, the spectral bands in ASTER GEDv3 first need to be spectrally

adjusted to the Landsat thermal bands. Landsat 4, 5 and 7 thermal sensor response (band 6) spans

~10.4-12.5μm, while that for Landsat 8 (band 10) spans ~10.5-11.5μm. The spectral adjustment

was accomplished using the ASTER narrow-band emissivity values from bands 13 (~10.6 μm)

and 14 (~11.3 μm) and building a regression function based on laboratory measurements (section

4.1).

Since the ASTER GEDv3 represents the average emissivity for 2000-2008, the emissivity

for any particular Landsat observation needs to be adjusted to capture changes due to vegetation

phenology. A methodology for this has been established by using Normalized Difference

Vegetation Index (NDVI) data derived from Landsat coupled with the emissivity and mean

NDVI products derived from ASTER GEDv3 to estimate surface emissivity at a particular

Landsat overpass time. A similar procedure is used to adjust for changes in snow cover using the

Normalized Difference Snow Index (NDSI) [58] derived from Landsat. The emissivity

10
adjustment does not account for changes in soil moisture (e.g. after rainfall), or changes in

ephemeral water.

4.1 Emissivity Spectral Adjustment

In order to use ASTER GEDv3 emissivities for the correction of Landsat thermal band(s),

they first need to be spectrally adjusted for Landsat thermal bands – band 6 (10.4-12.5 μm) for

Landsat 4 to 7, and band 10 for Landsat 8– using ASTER emissivities from bands 13 (10.6 μm)

and 14 (11.3 μm).

Figure 1 shows the spectral response functions of ASTER and Landsat thermal bands.

While there is no overlapping ASTER response for the longwave part of the Landsat 5/7 band

(11.75-12.5 µm), there is very low spectral contrast across this wavelength region, whereas there

is significant variation between 10 and 11.75 µm, allowing us to use a combination of the two

ASTER longwave bands 13 and 14. We estimated Landsat spectral emissivity associated with

band 6 of Landsat 4 to 7 using the broadband regression approach developed by Ogawa [59]:

𝜖𝜖10.4−12.5 = 𝑐𝑐13 𝜖𝜖13 + 𝑐𝑐14 𝜖𝜖14 + 𝑐𝑐 (4)

where c13, c14 and c are the regression coefficients, and 𝜖𝜖13 and 𝜖𝜖14 are ASTER GEDv3

emissivity values for bands 13 and 14, respectively.

The regression coefficients c13 and c14 in Eq. (4) are obtained using a set of rocks, soils,

vegetation, water and ice emissivity spectra from the ASTER spectral library [60], the Moderate-

Resolution Imaging Spectroradiometer (MODIS) spectral library

(http://www.icess.ucsb.edu/modis/EMIS/html/em.html) and sand samples from 10 sand dune

sites used to validate ASTER GED [57]. A total of 150 spectra were used to estimate 𝜖𝜖10.4−12.5,

𝜖𝜖13, and 𝜖𝜖14 using Eq. (4). The set of regression coefficients derived for each Landsat band’s

spectral response function is shown in Table 2. The root mean square error (RMSE) between

11
estimated and calculated emissivity values were 0.00085, 0.0013, 0.0011, and 0.0009 for Landsat

4, 5, 7 and 8, respectively (i.e. negligible). Note how the weighting of each ASTER band

depends on the amount of overlap with the Landsat thermal band, e.g. Landsat 7 and 8 have more

overlap with ASTER band 13 resulting in a larger coefficient value for that band, while 𝑐𝑐14 has

larger contribution for both Landsat 5 and 7 compared to Landsat 8.

Figure 2 shows an example cutout of spectrally adjusted emissivity derived from the

regression for Landsat 7 using the ASTER GEDv3 product. This region over the southwest USA

exhibits a large spatial heterogeneity in emissivity due to a mix of different land cover types

including water, bare soils, sand dunes, croplands, semi-arid shrublands and forest. Surface

emissivities range from approximately 0.94 to 0.99, with lower values occurring over bare

regions consisting mostly of quartz sands such as the Algodones Dunes and Kelso sand dunes

(Figure 2). The highest emissivity values occur over water (Salton Sea), densely vegetated crops

(Imperial Valley), and maffic rocks over the Mojave Desert.

The emissivity spectra for the Algodones and Kelso sand dunes sites are shown in Figure

3. In order to measure the emissivity of these sites, sand samples were collected and their

hemispherical reflectance measured using a Nicolet 520 Fourier transform infrared spectrometer.

The hemispherical reflectance was then converted to emissivity using Kirchhoff’s law (𝜖𝜖 = 1 –

r). The Algodones dune site contains mostly quartz sands resulting in the classic quartz doublet

feature between 8-9 microns, while Kelso dunes have a mix of quartz, feldspar, and potassium

with traces of magnetite [57]. Since the mineralogy and the emissivity of these sites are highly

constant over time, they serve as good validation targets for the ASTER GEDv3 spectral

adjustments for Landsat. Table 3 shows a comparison between lab-derived emissivities at both

sites and the spectrally adjusted ASTER GEDv3 emissivities for Landsat 4, 5, 7 and Landsat 8.

12
The absolute error represents the difference between the measured emissivity convolved to the

Landsat response function, and the ASTER-adjusted emissivity for each of the Landsat sensors

for the two sites. The Lab convolved and ASTER-adjusted values agrees to within ~1%. The

slightly higher emissivity values for adjusted ASTER GEDv3 are most likely due to sparse desert

grasses and reeds on dunes within the ASTER footprints which increases the emissivity, as

reported by Hulley et al. [57].

4.2 Vegetation Adjustment

Surface emissivity is strongly influenced by land cover type and land use [61]. Data from

the ASTER GEDv3 are able to resolve the surface emissivity of agricultural landscapes at field

scales (~100m), which is consistent with Landsat spatial resolutions. However, since the ASTER

GEDv3 represents average emissivity from 2000-2008, an emissivity adjustment is required to

modify this mean emissivity to land surface conditions at the time of the satellite overpass. In

the absence of soil moisture [57] emissivity over desert regions is largely invariant over time but

will increase with soil moisture [62], and vegetation cover [63]. Currently we do not account for

changes in soil moisture due to a lack of high resolution soil moisture data at field scale levels.

Several authors have proposed emissivity-NDVI relationships [24, 56, 64] based on the fact that

pixels with low NDVI values (bare soil, rocks, for example) usually have low spectral emissivity

values (<0.9), and pixels with high NDVI values, characteristic of photosynthetically active

vegetation, will have higher emissivity values (close to 1 for dense canopies). Therefore,

emissivity adjustments are necessary to estimate the emissivity changes due to annual and inter-

annual changes in plant phenology and vegetation density. This adjustment is usually relatively

smooth and gradual for natural vegetation covers, but can be sudden over agricultural areas due

to planting and harvesting activities.

13
 We used NDVI data derived from Landsat coupled with the emissivity and mean NDVI

products derived from ASTER GEDv3 to estimate the surface emissivity change at a particular

Landsat overpass time. The largest adjustments would occur for agricultural fields that were bare

for most ASTER observations, but were vegetated at Landsat overpass time, and vice-versa [65,

66].

First the bare soil component of the ASTER emissivity is estimated for each pixel by

apportioning the ASTER emissivity with the fractional vegetation cover estimated from the mean

NDVI product in the ASTER GEDv3 (Eq. 5). The emissivity is then adjusted based on the

Landsat derived fractional cover of the current observation (Eq. 5), adapted from Hulley et al.

[27].

𝜀𝜀𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿 = 𝑓𝑓𝑣𝑣,𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿 𝜀𝜀𝑣𝑣𝑣𝑣𝑣𝑣 + �1 − 𝑓𝑓𝑣𝑣,𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿 � 𝜀𝜀𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 (5)

𝜀𝜀𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴 − 𝜀𝜀𝑣𝑣𝑣𝑣𝑣𝑣 𝑓𝑓𝑣𝑣,𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴

where, 𝜀𝜀𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 =
1 − 𝑓𝑓𝑣𝑣,𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴

𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑚𝑚𝑚𝑚𝑚𝑚 − 𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁
where, 𝑓𝑓𝑣𝑣,𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿/𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴 = 1 −
𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑚𝑚𝑚𝑚𝑚𝑚 − 𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑚𝑚𝑚𝑚𝑚𝑚

where εASTER is the ASTER GEDv3 emissivity which has been spectrally adjusted using Eq. 4,

εLandsat is the computed Landsat emissivity, εbare is the estimated ASTER bare soil emissivity, εveg

is the vegetation emissivity, which is set to 0.98 (typical spectral library value for Landsat and

ASTER), fv,Landsat/ASTER is the vegetation cover fraction corresponding to Landsat or ASTER

NDVI as indicated by the superscript. The same methodology is used to adjust for changes in

snow cover by replacing NDVI with NDSI data. The total emissivity uncertainty (𝛿𝛿𝛿𝛿) associated

with this procedure is estimated as a combination of the standard deviation associated with the

ASTER GEDv3 emissivity product, 𝛿𝛿𝛿𝛿𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 (standard deviation of emissivity observations from

2000-2009, which inherently includes the uncertainty due to vegetation phenology), the RMSE

14
of the regression spectral fitting in Eq. 4, 𝛿𝛿𝛿𝛿𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , and the uncertainty associated with the TES

algorithm, 𝛿𝛿𝛿𝛿𝑟𝑟𝑟𝑟𝑟𝑟 , which was estimated as 0.0164 for band 13 and 0.0174 for band 14 [25]:

2 2 2
(𝛿𝛿𝛿𝛿𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 + 𝛿𝛿𝛿𝛿𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 + 𝛿𝛿𝜀𝜀𝑟𝑟𝑟𝑟𝑟𝑟 )
𝛿𝛿𝛿𝛿 = � (6)
3

The emissivity uncertainty is combined with an uncertainty model that estimates the uncertainty

associated with the atmospheric correction procedure detailed in Cook et al. [42]

Figure 4 shows an example of this methodology for a Landsat 7 cutout over agricultural

fields near Copeland, Kansas, USA. The NDVI values for the field are shown in Figure 4 (a &

b). According to the Kansas land use database (http://kars.ku.edu/maps/klcp2005/), the land

cover consists mostly of irrigated cornfields. Therefore, we chose dates illustrating the maximum

changes in emissivity from bare/fallow fields to fully matured crops during the summertime. The

usual planting date for irrigated corn in the area is April-May and the crops mature around July-

August before being harvested in September-October [67]. Figure 4 (a) shows fields mostly

consisting of bare soils with low NDVI values on April 24th 2002 before planting. Fixed-pivot

crop fields with high NDVI in 4 (a) represent matured winter wheat fields. Figure 4 (b) shows

increases in NDVI on Aug 14th when the corn starts to fully mature according to the USDA

planting and harvesting calendar [67]. Figure 4 (c) shows prior emissivity as in the baseline

ASTER GEDv3 since NDVI is mostly near zero across the cutout (bare, or senesced ground),

and 4 (d) final emissivity after fractional vegetation change based on Landsat NDVI. Figure 4 (e)

shows the percent emissivity change for Landsat 7 TIR band 6 using the methodology described

above in Eq. 5. Fully matured plants have high NDVI values >0.8 resulting in higher emissivity

values, which can be seen in Figure 4 (f) as an increase in baseline ASTER GEDv3 emissivity

values by as much as 3%. The 3% emissivity difference is equivalent to a temperature difference

15
of 2.2 K for a material at 300K, when measured by a sensor with filter centered around 11.3 µm

such as Landsat 7. Therefore, the vegetation emissivity adjustment associated with the algorithm

is an important step in ensuring accurate LST retrievals over landscapes with pronounced

temporal changes in vegetation cover.

4.3 Uncertainty Analysis

A full understanding of LST uncertainty and its sources is a critical aspect for meaningful

application and analysis of LST data over time. For example, an accurate knowledge of sources

of uncertainty is crucial for the analysis of long time series spanning several decades such as the

suite of Landsat sensors. There are several sources of uncertainty that are involved in the LST

retrieval e.g. uncertainty in the atmospheric profile, emissivity, water vapor or cloud proximity.

Uncertainty estimates can be a useful metric for users to investigate or mask pixels close to

clouds. Laraby et al. [68] noted that Landsat LST values were underestimated for pixels near

clouds, and reported progress towards developing an uncertainty estimate on LST retrievals due

to cloud proximity. The uncertainty associated with the atmospheric correction of Landsat

thermal data is described in detail in Cook et al. [42, 69]. In this section, we focus on emissivity

uncertainty and its effect on the LST retrieval for Landsat.

The effect of emissivity uncertainty on Landsat LST retrievals was estimated using the

Temperature Emissivity Uncertainty Simulator (TEUSim) developed by Hulley et al. [25].

TEUSim was developed to quantify the effects of algorithmic, atmospheric, and measurement

uncertainties on the retrieval of LST and emissivity from a number of sensors including MODIS,

Visible Infrared Imaging Radiometer Suite (VIIRS), ASTER, and Landsat TM/ETM+/TIR data.

The uncertainties are estimated with radiative transfer simulations (MODTRAN 5.2) using 382

16
global radiosonde profiles [70] and 155 surface emissivity spectra [60] as input for a range of

viewing angle configurations and surface temperatures for each profile.

To estimate the uncertainty associated with the ASTER GEDv3 emissivity correction for

Landsat, we assumed a perfect atmospheric correction (no assumed error in atmospheric

profiles), and used laboratory spectra of 26 rocks, 20 soils, 9 sands, and 4 graybodies (water, ice,

snow, conifer) from the ASTER spectral library with random uncertainty of 1.5% assigned to

each endmember in the forward simulation. This is consistent with the average uncertainty in the

ASTER GEDv3 product [27]. In addition, we assigned a standard fixed ‘classification’

emissivity for each class of spectra by assigning each class with the average emissivity of that

class: graybody = 0.9897, sands = 0.9733, soils = 0.9889, rocks = 0.9745. This is the more

traditional way of applying an emissivity correction per land cover type for single-channel

inversion method [24, 55, 56].

Figure 5 shows LST uncertainty distributions derived from TEUSim plotted versus Total

Column Water (TCW) and simulated LST for the Landsat LST single-channel inversion

algorithm using the ASTER GEDv3 (left) and emissivity classification (right). The LST

uncertainties using the physically retrieved ASTER GEDv3 are on average by a factor of 2

smaller than those of the classification approach i.e. with an RMSE of 1.23 K using ASTER

GEDv3 compared to 2.51 K for the classification approach. Note the cold bias of more than 2 K

for the classification approach, which is usually due to an overestimation of emissivity for each

class. Both approaches have the largest uncertainties under dry and warm conditions when the

surface has a much larger contribution to the observed radiance than the atmosphere, thereby

increasing the sensitivity to emissivity error.

17
5 LST Retrieval
Retrieval of LST from the Landsat series of TIR sensors using a single band approach

involves first calculating the surface emitted radiance using the estimated atmospheric

parameters (section 3), emissivity (section 4) and the observed Landsat radiance for band 6, and

then inverting the Planck function. First, Eq. (1) can be rearranged to calculate the surface

emitted radiance for the relevant Landsat sensor TIR band:

𝐿𝐿 − 𝐿𝐿↑
− (1 − 𝜖𝜖 ) 𝐿𝐿↓
𝐿𝐿𝑠𝑠 = 𝐵𝐵(𝑇𝑇𝑠𝑠 ) = 𝜏𝜏 (7)
𝜖𝜖

where 𝐵𝐵(𝑇𝑇𝑠𝑠 ) is the Planck function at temperature 𝑇𝑇𝑠𝑠 :

𝑐𝑐1 𝜆𝜆−5
𝐵𝐵(𝑇𝑇𝑠𝑠 ) = 𝑐𝑐2 (8)
�𝑒𝑒 𝜆𝜆𝑇𝑇𝑠𝑠 − 1�

where 𝑐𝑐1 = 1.19 × 10-16 W m-2 sr-1, and 𝑐𝑐2 = 1.44 × 104 µm K.

Theoretically the temperature can then be retrieved by inverting the Planck function with the

emitted radiance from Eq. (7) as follows:

𝑐𝑐2
𝑇𝑇𝑠𝑠 =
𝑐𝑐1 𝜆𝜆𝑐𝑐 −5 (9)
𝜆𝜆𝑐𝑐 ln �1 + 𝐿𝐿𝑠𝑠 �

where 𝑇𝑇𝑠𝑠 is the retrieved LST, and 𝜆𝜆𝑐𝑐 is the sensor’s central wavelength equivalent to a delta

function response. However, this formulation will increasingly become inaccurate for a sensor’s

spectral response that deviates from delta function behavior, for example the Landsat thermal

response functions covering a range of wavelengths in the TIR between 10—13µm (Figure 1).

Instead, we use a look up table (LUT) approach similar to the procedure used to calculate

ASTER brightness temperatures [71]. In this approach the Planck function is used to compute

18
expected radiances for each respective Landsat sensors’ spectral response over a range of

temperatures in 0.01 K intervals that encompass the full range of expected Earth-like

temperatures (typically 150 to 380 K). This results in a table of values of radiances versus

temperatures for each given sensor. The table can then simply be ‘inverted’ by interpolating to

get the retrieved temperature (LST), given the estimated surface emitted radiance (from Eq. 6).

The table can be modified to any desired precision by decreasing the step size interval (e.g. from

0.01 to 0.001 K).

Figure 6 shows an example of the LST retrieved using this approach for a Landsat 5

scene on 2 February, 2005 over the southwestern USA in the vicinity of the Salton Sea, CA. The

scene includes a wide range of surface types and temperatures such as water, managed

agricultural, natural vegetation and barren surfaces. The LST has a wide range from near

freezing over the peaks of the Anza Borrego’s in the west and southwest to more than 30ºC over

semi-arid shrublands surrounding the Salton Sea. It is important to note that although the thermal

bands are retrieved at lower resolution, e.g. 60 m for Landsat 7, the thermal bands are resampled

to 30 m pixel size in order to be consistent with the VSWIR multi-spectral bands [72].

6 LST Validation
6.1 SURFRAD sites

Surface Radiation Budget Network (SURFRAD) stations provide ground-based

observations of components of surface radiation budget e.g. upwelling and downwelling, solar

and infrared radiation etc. [73, 74] and have been used in the past to validate LST products from

ASTER and MODIS [75], and VIIRS [76]. The upwelling and downwelling longwave radiances

are measured by a downward and upward looking pyrgeometer, respectively. The uncertainty

associated with the field pyrgeometer is about 5 W m-2 [77], and an effective footprint of the

19
upwelling observation, deployed on a 10-m tower, is about 70x70 m2 [78]. The quality flagged

and aggregated daily SURFRAD data is available from an ftp site

(ftp://aftp.cmdl.noaa.gov/data/radiation/surfrad/). Data are available in 3 minutes interval before

Jan 2009, and every minute thereafter. All instruments are replaced on an annual basis with

freshly calibrated instruments. SURFRAD instruments are calibrated using three standards

maintained at NOAA's Field Test and Calibration Facility at Boulder, CO, which are traceable to

the World Radiometric Reference (WRR) at the World Radiation Center in Davos, Switzerland

[73]. There are seven SURFRAD stations, distributed across the United States in different

climatological regions. However, we excluded three of the SURFRAD sites based on surface

heterogeneity issues resulting in a mismatch between ground and satellite footprints and

described in more detail in the Appendix. Table 4 shows details of the four selected SURFRAD

sites for the present study.

The SURFRAD instrumentation measures the sum of reflected downwelling radiance,

and the surface leaving radiance, an integral of broad range of wavelengths (3.5-50 µm), dictated

by the Stefan-Boltzmann law.

(10)
𝐹𝐹 ↑ = (1 − 𝜖𝜖) 𝐹𝐹 ↓ + 𝜖𝜖𝑏𝑏 𝜎𝜎 𝐿𝐿𝐿𝐿𝐿𝐿 4

where, 𝐹𝐹 ↑ , 𝐹𝐹 ↓ are the measured upwelling, or downwelling longwave flux, respectively, 𝜖𝜖𝑏𝑏 is the

broadband longwave surface emissivity, 𝜎𝜎 is the Stefan-Boltzmann constant, and LST is the

surface skin temperature.

The LST is then found by inverting the surface-leaving component (Eq. 10):

0.25
1 ↑ ↓ (11)
𝐿𝐿𝐿𝐿𝐿𝐿 = � �𝐹𝐹 − (1 − 𝜖𝜖) 𝐹𝐹 ��
𝜖𝜖𝑏𝑏 𝜎𝜎

20
 The majority of the radiation emitted by terrestrial objects fall within the longwave

infrared (8-14 m) domain [79]. Therefore, the broadband emissivity (𝜖𝜖𝑏𝑏 ) is estimated from a

spectral-to-broadband regression relationship following Ogawa et al [80] using emissivities from

the ASTER GEDv3 product [27]:

(12)
𝜖𝜖𝑏𝑏 = 0.128 + 0.014 𝜖𝜖a10 + 0.145 𝜖𝜖a11 + 0.241 𝜖𝜖a12 + 0.467 𝜖𝜖a13 + 0.004 𝜖𝜖a14

where 𝜖𝜖a𝑗𝑗 , 𝑗𝑗 = 10 to 14 denote narrowband ASTER emissivities centered on 8.3 µm, 8.6

µm, 9.1 µm, 10.6 µm, and 11.3 µm, respectively. The SURFRAD measurements of longwave

radiation have been used in the past to validate LST products for ASTER and MODIS [75], and

VIIRS [78], for example.

For the validation of the Landsat LST product, we processed Landsat scenes between

2000 and 2012 (2003-2012 for Sioux Falls) for the respective path and row as listed in Table 4.

We excluded Landsat scenes with cloud fractions greater than 70% of the scene, as indicated in

the metadata (MTL) text file. The Landsat scenes are in UTM projection and are converted to

Latitude and Longitude grids in order to collocate with the SURFRAD station locations. Since,

the MTL text for each Landsat scene assigns the UTM Zone corresponding to the upper-left

corner, proper care must be taken to match the SURFRAD locations with the correct UTM Zone.

For example, GWN scenes with path 23, row 36 spans two UTM Zones 15 and 16. However,

these scenes are assigned to UTM Zone 15, while GWN site is located at the UTM Zone 16. This

was resolved by projecting both, the geo co-ordinates of SURFRAD site and the scene pixels, to

the same UTM Zone.

The Landsat LST product includes cloud mask information (cfmask) on a per pixel basis

[81]. Cloud contamination was minimized by removing data in which at least one pixel within a

5×5 window around the SURFRAD site was classified as cloud. Undetected clouds can result in

21
underestimation of retrieved Landsat LST by as much as 10 K. Any suspected cloud

contaminated matchups with SURFRAD were investigated qualitatively, and flagged for

removal. In order to reduce uncertainties due to spatial heterogeneity during daytime

observations, we removed observations that had a standard deviation exceeding 1 K for 3×3

pixels surrounding the site. All cloud-free processed data were then matched with the in-situ data

for the closest observation time.

In this study, we considered four of the most homogeneous SURFRAD stations at

Landsat ~100m pixel scales: Bondville (BON), Goodwin Creek (GWN), Penn State University

(PSU) and Sioux Falls (SXF) (Table 4). The results are summarized by the bias and the root

mean square error (RMSE) between Landsat LST and SURFRAD LST matchups.

6.1.1 Bondville station (BND)

The BON station (Figure 7e) is located in a grass field, surrounded by corn and soybean

fields that have typical planting dates from Apr-May. Figure 8(a) shows a scatterplot between

Landsat retrieved LST and the in-situ skin temperature calculated using Eq. 10. The bias

(RMSE) in LST are 1.3 K (2.5 K), and 1.1 K (2.5 K) for Landsat 5 and 7, respectively.

To better understand the nature of the bias for both sensors at these sites, we illustrate the bias

(Landsat minus SURFRAD LST) by month using box plots [ 1] in Figure 9. Figure 9(a) shows the

seasonal variability in the bias distribution between the retrieved LST and SURFRAD data.

During the growing season the bias is highest (~4 K, around May) when there is larger soil

exposure between newly planted crops and warmer effective temperatures over the larger ground

footprint of Landsat when compared to the SURFRAD sensor. The LST Standard deviation

1 The central mark in box plot indicates the median, and the bottom and top edges of the box indicate the 25th and
75th percentiles, respectively. The outliers are indicated using the ‘+’, or ‘o’ symbols for Landsat 5 and Landsat 7,
respectively.

22
(STD) of 3x3 Landsat pixels surrounding the site also has higher variability (>1K) during the

month of May for both L5 and L7. The bias decreases to a minimum during the summer as the

crops fully mature in Jul-Aug resulting in cooler and more homogeneous temperature

distribution around the site. The bias increases again during the harvesting season (Sept-Oct,

~2K) due to difference in the exposed soil within Landsat and SURFRAD footprints.

6.1.2 Goodwin Creek Station (GWN)

The GWN station is located in a rural pasture grassland surrounded by broadleaf

deciduous forest (Figure 7f). The site is largely homogeneous for a 100 m radius around the site

with few trees and mostly grasses. The nearest stand of trees around the site are ~60 m away, and

there is a watershed ~120 m to the west. A spatial variation analysis of a 3x3 window of the

Landsat LST around the site did not show significant seasonal variability and as a result the site

can be regarded as homogeneous at the Landsat pixel resolution size. Figure 8(b) shows that the

retrieval performed satisfactorily with bias (RMSE) 0.5 K (2.0 K) and 1.1 K (2.3 K) for Landsat

5 and 7, respectively.

Figure 9 (b) shows a larger bias and variance during Jul-Aug for both Landsat 5 and 7

retrievals are higher, and this can likely be attributed to summer monsoonal conditions when

atmospheric transmissivity decreased to less than 0.7 during Jun-Aug, resulting in larger

uncertainties in the atmospheric correction and LST single-channel retrieval.

6.1.3 Penn State Station (PSU)

The PSU SURFRAD station is located on an agricultural research farm in an area

surrounded by mostly corn and soybean fields (Figure 7c) , with similar planting dates starting in

May. Figure 8(c) shows a comparison of Landsat LST with the ground station data at the PSU

site with a bias (RMSE) of 0.2 K (2.3 K), and 0.3 K (2.1 K) for Landsat 5 and Landsat 7,

respectively. However, similar to the BND site, the LST differences have a seasonal nature and

23
increase to 3-4 K for both sensors during the planting season from May-July (Figure 9c) most

likely due to discrepancies and changes in land cover between the ground sensor field of view

and the satellite pixel resolution. This site is also located in an area with significant atmospheric

moisture during the summertime, which also results in higher uncertainties in atmospheric

correction and the LST single-band retrieval.

6.1.4 Sioux Falls Station (SXF)

The SXF station is located on a grassland about 130 m southeast of the Earth Resources

Observation and Science (EROS) Data Center, SD (Figure 7b). The site is homogeneous at the

Landsat resolutions of ~100m, but is not suitable for validation of moderate resolution LST (~1

km) due to the close proximity of surrounding crop fields and heterogeneous mixture of land

cover including buildings, roads, and other structures.

Figure 8(d) shows the comparison between Landsat retrieved LST and the ground station

data with a bias (RMSE) of 0.7 K (2.1 K), and 1.1 K (2.6 K) for Landsat 5 and Landsat 7,

respectively. Similar to previous sites, Figure 9(d) shows that higher biases are observed during

summer months for both sensors (2-4 K), while the winter months have fewer observations due

to persistent cloud cover.

6.1.5 SURFRAD Validation Summary

Figure 10 shows an average histogram distribution of the LST differences between

Landsat 5/7 and SURFRAD data. The overall bias (RMSE) is 0.7 K (2.2 K) and 0.9 K (2.3 K) for

Landsat 5 and Landsat 7, respectively. About 5-7 % of the data have a bias greater than +4 K.

These results indicate an accuracy within 1 K and total uncertainty of about 2.5 K RMSE for the

Landsat 5 and 7 LST algorithm over the SURFRAD land sites for a wide range of temperatures

and clear sky conditions.

24
 The positive bias can largely be attributed due to changes in phenology during the

growing and harvesting seasons at the agricultural sites (e.g. BND, and PSU in Figure 9) mostly

as a result of the different mixture of surface covers observed by Landsat and SURFRAD

sensors. This phenomenon was also observed and discussed by Guillevic et al. [78] when

validating VIIRS LST at these sites. The result is that the Landsat observations are more

sensitive to changes in phenology in the surrounding crop fields resulting in higher summertime

biases that are exacerbated by higher spatial temperature variability.

The positive biases were also higher during the summer months at grassland sites (e.g.

GWN, and SXF) most likely due to a combination of higher temperatures (driving spatial

variability) and humid conditions from summertime monsoonal moisture that increases

atmospheric correction uncertainty.

6.2 Inland Water Sites

We also validated the Landsat LST product over two inland water bodies in

California/Nevada (viz. Lake Tahoe and Salton Sea). Water has high and very well-known

emissivity (𝜖𝜖), and compared to land surfaces, temperatures are spatially homogeneous and

change slowly, making it a good calibration and validation target for thermal infrared data [82-

84]. Since the ASTER GED prescribes a laboratory-derived emissivity over inland water

surfaces using a land/water mask, the results here mostly describe the accuracy of the

atmospheric correction technique. We collocated Landsat overpass scenes over the two inland

water sites, and filtered them for clouds using the techniques described in previous section.

6.2.1 Lake Tahoe

Lake Tahoe is a large inland freshwater body located in the Sierra Nevada Mountains at

39°N and 120°W at the California-Nevada border. It has a surface elevation of 1898 m and is

25
about 35 km long and 19 km wide [84]. It is regarded as the sixth largest lake by volume, behind

the five Great Lakes, and is the second deepest in the United States after Crater Lake in Oregon.

JPL has been continuously measuring skin and bulk temperatures every 2 minute on four

automatic buoy sites since 1999 [85]. Each buoy consists of a custom-built radiometer developed

by JPL, which are calibrated every 10 minute, and are accurate to ±0.2 K [84]. The radiometer is

mounted on a pole, extended beyond the buoy, at about 1 m above the surface of water. The

radiometric measurements are converted to skin temperatures by accounting for the effects of

emissivity and reflected downwelling sky radiation. The emissivity of water is obtained from the

ASTER spectral library [60], while the downwelling radiance is computed using a radiative

transfer model (MODTRAN) with NCEP input data [84].

6.2.2 Salton Sea

The second water validation site, Salton Sea, is located in the southern California at 33.3°N, and

115.8°W, and is a shallow, saline lake with about 900 km2 surface area and surface elevation 71

m below sea level. The JPL radiometer is situated at the southwest corner of the lake on a

platform, and has been taking measurements since the end of 2007. The combination of two lake

sites provide contrasting atmospheric conditions with surface temperatures ranging from ~4 to

35°C, and wide variety of atmospheric conditions [86].

6.2.3 Lake Validation Results

For the validation, we collocated Landsat overpasses over each site, and compared the

satellite-derived LST with collocated and temporally matched to 2 minutes temperature retrievals

from in situ radiometers for the time range of 2000-2011 for Lake Tahoe, and 2007-2011 for

Salton Sea. Figure 11 shows the validation matchup results for clear sky Landsat LST data

retrieved over Salton Sea and Lake Tahoe. The mean bias (RMSE) between the retrieved

Landsat LST and in-situ skin temperature was − 0.2 K (0.6 K) and 0.4 K (0.8 K) for Landsat 5

26
and Landsat 7, respectively (Figure 11 (a)). Similarly, the mean bias (RMSE) for the Salton Sea

observations was −0.4 K (0.7 K) and 0.3 K (0.5 K) for Landsat 5 and Landsat 7, respectively

(Figure 11 (b)). Figure 12 shows an average histogram distribution of the difference between the

retrieved Landsat 5/7 and in situ data for both lakes. The results show overall good agreement of

<1 K with a mean bias (RMSE) for Landsat 5 and Landsat 7 of -0.3 K (0.6 K), and 0.4 K (0.7 K),

respectively.

7 Conclusions
Due to lack of an accurate retrieval retrieval methodology for land surface temperature

(LST) at global scales and a suitable surface emissivity dataset, up to now thermal infrared (TIR)

data from Landsat satellites have been underutilized. To address these issues, we a team from

JPL and RIT have developed an operational single-channel (SC) algorithm using the ASTER

Global Emissivity Dataset (GEDv3) that will be used to generate standardized LST product for

the entire Landsat series extending back to 1982.

The developed single channel LST retrieval methodology for Landsat first requires an

atmospheric correction of the TIR band, followed by an emissivity correction and a single-

channel (SC) inversion via a Look-Up-Table (LUT) approach. The atmospheric correction is

performed using the MODTRAN 5.2 radiative transfer model with numerical weather model

data from NARR (3-hourly, ~32 km) for North America, and MERRA-2 (6-hourly, ~50 km),

while emissivity correction uses the ASTER GEDv3 at 100m spatial resolution. The physically

retrieved emissivity values provided by ASTER GEDv3 enable a more accurate and consistent

emissivity correction of Landsat thermal data at its native resolution. Spectral adjustment of the

ASTER GEDv3 for Landsat thermal bands and dynamic adjustments accounting for vegetation

phenology and snow cover are discussed in detail including uncertainty analysis of the emissivity

27
effects on the retrieved LST. The spectral adjustments were validated for the Landsat suite of

sensors at two pseudo invariant sand dune sites with an average 0.54% absolute error. For a

Landsat 7 scene over corn fields near Copeland, Kansas, USA, we qualitatively demonstrated an

increase in emissivity due to phenology from baseline ASTER GEDv3 by up to 3%.

The SC retrieval technique was validated over land at four SURFRAD sites using

observations from 2000-2012. On average, the retrieved Landsat LST had an accuracy of ~1 K

and total uncertainty of 2.5 K. The monthly average distribution of LST error at each site (Fig. 9)

revealed that biases at the agricultural sites (BND, PSU) were higher (Landsat warmer than

ground station) during the spring/summer planting and harvesting seasons. This is explained by

the different surface areas observed by Landsat and the ground sensors, resulting in a Landsat

‘contamination’ from surrounding warmer exposed bare surfaces, whereas the ground sensors

observe grass. Other contributing factors are higher summertime temperatures and higher spatial

variability in surface temperatures, combined with higher water vapor loadings that result in

larger uncertainties in the atmospheric correction. The Landsat LST retrievals were also

validated over two inland waterbodies using collocated in-situ skin temperature data from Lake

Tahoe, CA/NV between years 2000-2011, and the Salton Sea, CA from 2008-2011. The results

showed very good agreement with in situ data with a bias of −0.3 K and an RMSE of 0.6 K for

Landsat 5, and a bias of 0.4 K and RMSE 0.7 K for Landsat 7. Since the emissivity of water is

well known, this validation exercise represents mostly an evaluation of the atmospheric

correction algorithm.

The close agreement (<1 K) between Landsat and in-situ data at the lake sites indicate that

the larger uncertainties over the land SURFRAD sites (2.5 K) are driven by the large spatial

heterogeneity of surface temperature and vegetation phenology changes within the Landsat field

28
of view, which are not represented by the ground observations. These results strongly support

previous studies [75] that the SURFRAD sites, while suitable for validation of moderate

resolution sensors (ASTER, Landsat) with careful consideration of sensor/site geolocation (see

discussion in Appendix), precise ground sensor location (via communication with site PI’s), and

time of year, are not suitable for validating moderate resolution thermal data from MODIS (1km)

and VIIRS (750m), unless validation is performed with nighttime only data to minimize spatial

variability, or a model is used to account for spatial heterogeneity [78]. A more extensive

validation on global scale over a wider variety of dedicated land sites will be presented in the

near future in order to reach Stage-2 validation status..

The operational Landsat LST algorithm has been implemented at EROS LPDAAC and

processing of the Landsat 5 (TM) LST product over the Continental US will begin during the

Fall 2017 and made available through Earth Explorer (http://earthexplorer.usgs.gov). Evaluation

of the global LST product using MERRA-2 data for the atmospheric correction is ongoing. The

operational algorithm will, for the first time, produce a long-term LST record from the Landsat

thermal data record dating back to Landsat-4 (1982-present) using a consistent retrieval

methodology and emissivity correction method. This will enable the monitoring of long-term

temperature trends, utilization in surface energy balance models to estimate important drought-

monitoring indicators such as evapotranspiration and evaporative stress indices, detection of land

cover land use changes, water resource management over agricultural systems, and studies

related to the urban heat island effect.

8 Appendix
Validation of LST data from space requires a site that is homogeneous in temperature at

the scale of the imagery, ideally allowing several image pixels to be validated over the target site.

29
Subpixel land cover heterogeneity is one of the largest sources of uncertainty in the validation of

LST using ground sensors, particularly for moderate scale resolution sensors (1 km). Since the

Landsat sensors validated in this work (L5, L7) have only one thermal band, multispectral

Radiance-based (R-based) validation techniques [87-89] were not possible, leaving only

Temperature-based (T-based) validation methods as the primary source. For validation of LST

over water sites (e.g. Tahoe, Salton Sea) subpixel heterogeneity is not an issue. However,

because the emissivity of water is known a priori, the validation only tells us something about

the accuracy of the atmospheric correction. For validation over land sites, the NOAA SURFRAD

sites have been used frequently in past validation studies to validate LST data from MODIS,

VIIRS, ASTER [75, 78, 90]. SURFRAD sites are attractive because they are well maintained,

have good accuracy, and data are collected continuously, which allow close temporal matchups

with satellite data. However, their biggest drawback is spatial heterogeneity in cover types (and

hence temperature) around the sites as previously reported by Wang and Liang [75] and

Guillevic et al. [91]. This appendix focuses on assessment of this variability at particular sites,

justification for not using these sites in the validation study, and guidance for future validation

efforts at these sites.

Landsat data in the visible-shortwave infrared (VSWIR) are generally acquired at higher

spatial resolution than the thermal infrared (TIR). This is because TIR wavelengths are ~10x

longer than in the visible, necessitating 10x larger optics to achieve the same ground resolution

as in the visible. Therefore, Landsat 7 VSWIR data are acquired at 30m, while the TIR data are

acquired at 60m, which requires the thermal bands to be resampled to 30 m pixel size to be

consistent with the VSWIR bands [72]. Assuming that the resampled pixels do not introduce

additional error, spatial misrepresentation between ground sensor effective footprint and Landsat

30
pixels is a major contributor to validation uncertainty. At the SURFRAD sites, this subpixel

variability needs to be carefully assessed, and could be due to different surface cover fraction

surrounding the site or different vegetation phenologies at the sites surrounded by

crop/grasslands. Based on our investigations of site variability, we decided to remove the

following sites from our study to reduce unwanted biases due to satellite/ground sensor footprint

mismatches and surface temperature variability: DRA, FPK, and TBL.

8.1.1 Desert Rock (DRA)

The DRA site (36.6232, -116.01962) is located in an arid valley-about 90 km northwest

of Las Vegas, NV (Figure 13). The site consists mostly of a mix between small creosote shrubs

and exposed gravel and bare soils, but also evident are patches of dark gravel (black areas) that

were dumped in the vicinity of the site at some point in the past (personal communication, Gary

Hodges and John Augustine, Site PI’s, 04/2013). As an example of subpixel variation, in Figure

13, we show retrieved LST from the Hyperspectral Thermal Emission Spectrometer (HyTES)

airborne sensor at ~6 m spatial resolution at the DRA site (available at https://hytes.jpl.nasa.gov,

ID: 2014-07-03 21:16:06 ). Details of the HyTES mission and algorithm are beyond the scope of

this paper. Interested readers are referred to [92]. Figure 13 (a) shows a Google Earth image of

the DRA site along with the footprints of L5/L7 TIR sensors compared with the ground sensor’s

field of view represented by a circle. Figure 13 (b) shows that the corresponding HyTES

retrieved LST over the darker material is up to 15 K warmer than over the surrounding areas.

Figure 13 (c), and (d) show the light colored and darker gravels, respectively. Since the L5 TIR

sensor footprint width is twice as large (120m) than that of the ground sensor (~60m), the

effective temperature would include the hotter darker material, which explains the persistent

positive bias at DRA of ~2.5 K for both L5 and L7 sensors during the daytime, and justifies our

31
decision for excluding this site. For other coarser sensors such as MODIS/VIIRS at kilometer

scale resolution, these biases may not be evident because of compensation from other areas such

as the cooler buildings and increased vegetation to the north (as shown in Figure 13b),

nevertheless that still does not justify the use of this site for coarser resolution sensors.

8.1.2 Fort Peck (FPK) station

The FPK site (48.307933, -105.101729) is located in a shrub-grassland in the Fort Peck

Tribes Reservation, inside a small fence-protected area (25 m x 25 m). We observed large

seasonal differences between Landsat and in-situ LST data (~5 K for both sensors) at the FPK

site particularly during the summertime. During winter, there were very few observations since

the site experienced persistent cloud cover. The large summertime biases may be driven by the

grazing patterns of bison herds, which could result in significant contrast in the land cover inside

and outside the fence as discovered in previous studies [93]. Due to these reasons we excluded

this site from our validation study.

8.1.3 Table Mountain (TBL)

The TBL site (40.12559, -105.237792) is located on the plateau of Table Mountain, about

4 km east of the Rocky Mountain foothills, and about 10 km north of Boulder, CO, in a grassland

consisting of sparse grass and shrubs with exposed sandy mix of exposed rocks. The site showed

a consistent positive bias of ~3.5 K for L5/L7 as reported by previous authors when validating

ASTER LST products [75]. One hypothesis is that the location of the instrument near the edge of

a steep plateau results in warmer temperatures due to warmer south facing slopes from solar

insolation in the daytime hours at the satellite effective pixel scale. Therefore, for this study, we

did not include the SURFRAD site located at Table Mountain, Boulder Colorado.

32
 Acknowledgements. The research described in this paper was carried out at the Jet

Propulsion Laboratory, California Institute of Technology, government sponsorship

acknowledged. NKM would like to thank Nathan Healey and Kerry Cawse-Nicholson for their

comments on draft version of the paper. We would also like to thank SURFRAD PIs for

maintaining the stations, and in particular John Augustine and Gary Hodges for providing

images, and samples from the Desert Rock site. The ASTER GEDv3 and GEDv4 are both

available for download at

https://lpdaac.usgs.gov/dataset_discovery/community/community_products_table.

9 References
[1] J. Lloyd and J. Taylor, "On the temperature dependence of soil respiration," Functional
ecology, pp. 315-323, 1994.
[2] B. Bond-Lamberty and A. Thomson, "Temperature-associated increases in the global soil
respiration record," Nature, vol. 464, no. 7288, pp. 579-582, 2010.
[3] J. D. Kalma, T. R. McVicar, and M. F. McCabe, "Estimating land surface evaporation: A
review of methods using remotely sensed surface temperature data," Surveys in
Geophysics, vol. 29, no. 4-5, pp. 421-469, 2008.
[4] Z. Wan, P. Wang, and X. Li, "Using MODIS land surface temperature and normalized
difference vegetation index products for monitoring drought in the southern Great Plains,
USA," International Journal of Remote Sensing, vol. 25, no. 1, pp. 61-72, 2004.
[5] M. C. Anderson, C. R. Hain, B. Wardlow, J. R. Mecikalski, and W. P. Kustas,
"Evaluation of a drought index based on thermal remote sensing of evapotranspiration
over the continental U.S.," Journal of Climate, vol. 24, pp. 2025-2044, 2011.
[6] M. C. Anderson, J. M. Norman, J. R. Mecikalski, J. A. Otkin, and W. P. Kustas, "A
climatological study of evapotranspiration and moisture stress across the continental
United States based on thermal remote sensing: 2. Surface moisture climatology," (in
English), Journal of Geophysical Research-Atmospheres, vol. 112, no. D11112, Jun 9
2007.
[7] R. G. Allen et al., "Satellite-based energy balance for mapping evapotranspiration with
internalized calibration (METRIC) - Applications," (in English), Journal of Irrigation
and Drainage Engineering-Asce, vol. 133, no. 4, pp. 395-406, 2007/07//undefined 2007.
[8] O. Taconet, R. Bernard, and D. Vidalmadjar, "Evapotranspiration over an Agricultural
Region Using a Surface Flux Temperature Model Based on Noaa-Avhrr Data," (in
33
 English), Journal of Climate and Applied Meteorology, vol. 25, no. 3, pp. 284-307, Mar
1986.
[9] S. L. Harlan, A. J. Brazel, L. Prashad, W. L. Stefanov, and L. Larsen, "Neighborhood
microclimates and vulnerability to heat stress," Social Science & Medicine, vol. 63, no.
11, pp. 2847-2863, 12// 2006.
[10] B. Dousset et al., "Satellite monitoring of summer heat waves in the Paris metropolitan
area," International Journal of Climatology, vol. 31, no. 2, pp. 313-323, 2011.
[11] S. Flasse, S. Connor, A. Perryman, M. Thomson, and W. H. Organization, "The
contribution of satellite derived information to malaria stratification, monitoring and
early warning," 1997.
[12] M. O. Ruiz et al., "Local impact of temperature and precipitation on West Nile virus
infection in Culex species mosquitoes in northeast Illinois, USA," Parasit Vectors, vol. 3,
no. 1, p. 19, 2010.
[13] K. P. Willett, "Progress Report for the International Surface Temperature Initiative," pp.
ISTI Steering Committee, http://www.wmo.int/pages/prog/gcos/SC-
XIX/12%20Intl%20Surface%20Temp%20Initiative--Willett.pdf, 2011.
[14] NASA, "Exploring our Planet for the Benefit of Society, NASA Earth Science and
Applications from Space, Strategic
Roadmap,http://images.spaceref.com/news/2005/earth_roadmap.pdf," 2005.
[15] GCOS, "The Second Report on the Adequacy of the Global Observing Systems for
Climate in Support of the UNFCCC," pp. GCOS-82, April 2003 (WMO/TD No. 1143),
2003.
[16] R. De Fries, M. Hansen, J. Townshend, and R. Sohlberg, "Global land cover
classifications at 8 km spatial resolution: the use of training data derived from Landsat
imagery in decision tree classifiers," International Journal of Remote Sensing, vol. 19,
no. 16, pp. 3141-3168, 1998.
[17] J. A. Barsi, J. R. Schott, F. D. Palluconi, and S. J. Hook, "Validation of a web-based
atmospheric correction tool for single thermal band instruments," in Optics & Photonics
2005, 2005, pp. 58820E-58820E-7: International Society for Optics and Photonics.
[18] J. A. Barsi, J. R. Schott, S. J. Hook, N. G. Raqueno, B. L. Markham, and R. G.
Radocinski, "Landsat-8 Thermal Infrared Sensor (TIRS) vicarious radiometric
calibration," Remote Sensing, vol. 6, no. 11, pp. 11607-11626, 2014.
[19] USGS. (2014, October 7, 2016). Landsat 8 (L8) Operational Land Imager (OLI) and
Thermal Infrared Sensor (TIRS).
[20] D. P. Roy et al., "Landsat-8: Science and product vision for terrestrial global change
research," Remote Sensing of Environment, vol. 145, pp. 154-172, 2014.
[21] J. A. Sobrino, J. C. Jiménez-Muñoz, and L. Paolini, "Land surface temperature retrieval
from LANDSAT TM 5," Remote Sensing of Environment, vol. 90, no. 4, pp. 434-440,
4/30/ 2004.
[22] J. C. Jimenez-Munoz, J. Cristobal, J. A. Sobrino, G. Soria, M. Ninyerola, and X. Pons,
"Revision of the Single-Channel Algorithm for Land Surface Temperature Retrieval
From Landsat Thermal-Infrared Data," (in English), Ieee Transactions on Geoscience
and Remote Sensing, Article vol. 47, no. 1, pp. 339-349, Jan 2009.
[23] J. C. Jimenez-Munoz and J. A. Sobrino, "A Single-Channel Algorithm for Land-Surface
Temperature Retrieval From ASTER Data," (in English), Ieee Geoscience and Remote
Sensing Letters, vol. 7, no. 1, pp. 176-179, Jan 2010.

34
[24] W. C. Snyder, Z. Wan, Y. Zhang, and Y. Z. Feng, "Classification-based emissivity for
land surface temperature measurement from space," (in English), International Journal of
Remote Sensing, vol. 19, no. 14, pp. 2753-2774, Sep 20 1998.
[25] G. C. Hulley, C. G. Hughes, and S. J. Hook, "Quantifying uncertainties in land surface
temperature and emissivity retrievals from ASTER and MODIS thermal infrared data,"
(in English), Journal of Geophysical Research-Atmospheres, vol. 117, Dec 12 2012.
[26] J. A. Barsi, J. L. Barker, and J. R. Schott, "An atmospheric correction parameter
calculator for a single thermal band earth-sensing instrument," in Geoscience and Remote
Sensing Symposium, 2003. IGARSS'03. Proceedings. 2003 IEEE International, 2003, vol.
5, pp. 3014-3016: IEEE.
[27] G. C. Hulley, S. J. Hook, E. Abbott, N. Malakar, T. Islam, and M. Abrams, "The ASTER
Global Emissivity Dataset (ASTER GED): Mapping Earth's emissivity at 100 meter
spatial scale," (in English), Geophysical Research Letters, vol. 42, no. 19, pp. 7966-7976,
Oct 16 2015.
[28] G. C. Hulley and S. J. Hook, "The North American ASTER Land Surface Emissivity
Database (NAALSED) Version 2.0," (in English), Remote Sensing of Environment,
Article vol. 113, no. 9, pp. 1967-1975, Sep 2009.
[29] Z. M. Wan and Z. L. Li, "A physics-based algorithm for retrieving land-surface
emissivity and temperature from EOS/MODIS data," (in English), Ieee Transactions on
Geoscience and Remote Sensing, vol. 35, no. 4, pp. 980-996, Jul 1997.
[30] S. J. Hook, A. R. Gabell, A. A. Green, and P. S. Kealy, "A Comparison of Techniques for
Extracting Emissivity Information from Thermal Infrared Data for Geologic Studies," (in
English), Remote Sensing of Environment, vol. 42, no. 2, pp. 123-135, Nov 1992.
[31] P. S. Kealy and S. Hook, "Separating temperature & emissivity in thermal infrared
multispectral scanner data: Implication for recovering land surface temperatures.," IEEE
Transactions on Geoscience and Remote Sensing, vol. 31, no. 6, pp. 1155-1164, 1993.
[32] J. C. Jimenez-Munoz, J. A. Sobrino, C. Mattar, G. Hulley, and F. M. Gottsche,
"Temperature and Emissivity Separation From MSG/SEVIRI Data," (in English), Ieee
Transactions on Geoscience and Remote Sensing, vol. 52, no. 9, pp. 5937-5951, Sep
2014.
[33] T. Schmugge, A. French, J. C. Ritchie, A. Rango, and H. Pelgrum, "Temperature and
emissivity separation from multispectral thermal infrared observations," (in English),
Remote Sensing of Environment, vol. 79, no. 2-3, pp. 189-198, 2002/02//undefined 2002.
[34] A. Gillespie, S. Rokugawa, T. Matsunaga, J. S. Cothern, S. Hook, and A. B. Kahle, "A
temperature and emissivity separation algorithm for Advanced Spaceborne Thermal
Emission and Reflection Radiometer (ASTER) images," (in English), IEEE Transactions
on Geoscience and Remote Sensing, vol. 36, no. 4, pp. 1113-1126, Jul 1998.
[35] J. C. Jimenez-Munoz and J. A. Sobrino, "A generalized single-channel method for
retrieving land surface temperature from remote sensing data," (in English), Journal of
Geophysical Research-Atmospheres, Article vol. 108, no. D22, p. 9, Nov 2003, Art. no.
4688.
[36] Z. L. Li et al., "Satellite-derived land surface temperature: Current status and
perspectives," (in English), Remote Sensing of Environment, Review vol. 131, pp. 14-37,
Apr 2013.

35
[37] Z. M. Wan and J. Dozier, "A generalized split-window algorithm for retrieving land-
surface temperature from space," (in English), IEEE Transactions on Geoscience and
Remote Sensing, vol. 34, no. 4, pp. 892-905, Jul 1996.
[38] C. Coll and V. Caselles, "A split-window algorithm for land surface temperature from
advanced very high resolution radiometer data: Validation and algorithm comparison,"
(in English), Journal of Geophysical Research-Atmospheres, vol. 102, no. D14, pp.
16697-16713, Jul 27 1997.
[39] J. C. Price, "Land surface temperature measurements from the split window channels of
the NOAA 7 Advanced Very High Resolution Radiometer," Journal of Geophysical
Research, vol. 89, no. D5, pp. 7231-7237, 1984.
[40] Z. Qin, A. Karnieli, and P. Berliner, "A mono-window algorithm for retrieving land
surface temperature from Landsat TM data and its application to the Israel-Egypt border
region," (in English), International Journal of Remote Sensing, vol. 22, no. 18, pp. 3719-
3746, Dec 2001.
[41] F. Q. Li et al., "Deriving land surface temperature from Landsat 5 and 7 during
SMEX02/SMACEX," (in English), Remote Sensing of Environment, vol. 92, no. 4, pp.
521-534, Sep 30 2004.
[42] M. Cook, "Atmospheric Compensation for a Landsat Land Surface Temperature
Product," Ph.D., Chester F. Carlson Center for Imaging Science, Rochester Institute of
Technology, Rochester, 2014.
[43] F. X. Kneizys et al., "The MODTRAN 2/3 Report & LOWTRAN 7 Model, F19628-91-
C-0132," presented at the Phillips Lab, Hanscom AFB, MA, 1996.
[44] A. Berk, "MODTRAN: A moderate resolution model for LOWTRAN 7," Spectral
Sciences, Inc., Burlington, MA., 1989.
[45] A. Berk et al., "MODTRANTM 5, A Reformulated Atmospheric Band Model with
Auxiliary Species and Practical Multiple Scattering Options: Update," presented at the in
Proc SPIE, Algorithms and Technologies for Multispectral, Hyperspectral, and
Ultraspectral Imagery XI, Bellingham, WA, 2005.
[46] J. A. Barsi et al., "Landsat TM and ETM+ thermal band calibration," Canadian Journal
of Remote Sensing, vol. 29, no. 2, pp. 141-153, 2003.
[47] F. Mesinger et al., "North American regional reanalysis," (in English), Bulletin of the
American Meteorological Society, vol. 87, no. 3, pp. 343-+, Mar 2006.
[48] E. J. Becker, E. H. Berbery, and R. W. Higgins, "Understanding the characteristics of
daily precipitation over the United States using the North American Regional
Reanalysis," Journal of Climate, vol. 22, no. 23, pp. 6268-6286, 2009.
[49] M. S. Bukovsky and D. J. Karoly, "A brief evaluation of precipitation from the North
American Regional Reanalysis," (in English), Journal of Hydrometeorology, Article vol.
8, no. 4, pp. 837-846, Aug 2007.
[50] A. Molod, L. Takacs, M. Suarez, and J. Bacmeister, "Development of the GEOS-5
atmospheric general circulation model: evolution from MERRA to MERRA2," Geosci.
Model Dev., vol. 8, no. 5, pp. 1339-1356, 2015.
[51] M. G. Bosilovich, J. Y. Chen, F. R. Robertson, and R. F. Adler, "Evaluation of global
precipitation in reanalyses," (in English), Journal of Applied Meteorology and
Climatology, vol. 47, no. 9, pp. 2279-2299, Sep 2008.
[52] M. M. Rienecker et al., "MERRA: NASA's modern-era retrospective analysis for
research and applications," Journal of Climate, vol. 24, no. 14, pp. 3624-3648, 2011.

36
[53] E. DAAC, "GTOPO 30 Database. EROS Data Center Distributed Active Archive
Center," US Geological Survey, EROS Data Center, Sioux Falls, South Dakota, 1996.
[54] D. Shepard, "A two-dimensional interpolation function for irregularly-spaced data," in
Proceedings of the 1968 23rd ACM national conference, 1968, pp. 517-524: ACM.
[55] J. A. Sobrino et al., "Land surface emissivity retrieval from different VNIR and TIR
sensors," IEEE Transactions on Geoscience and Remote Sensing, vol. 46, no. 2, pp. 316-
327, 2008.
[56] E. Valor and V. Caselles, "Mapping land surface emissivity from NDVI: Application to
European, African, and South American areas," (in English), Remote Sensing of
Environment, vol. 57, no. 3, pp. 167-184, Sep 1996.
[57] G. C. Hulley, S. J. Hook, and A. M. Baldridge, "Validation of the North American
ASTER Land Surface Emissivity Database (NAALSED) Version 2.0 using Pseudo-
Invariant Sand Dune Sites," Remote Sensing of Environment, vol. 113, pp. 2224-2233,
2009.
[58] D. K. Hall, G. A. Riggs, and V. V. Salomonson, "Development of Methods for Mapping
Global Snow Cover Using Moderate Resolution Imaging Spectroradiometer Data," (in
English), Remote Sensing of Environment, vol. 54, no. 2, pp. 127-140, Nov 1995.
[59] K. Ogawa, "Mapping Surface Broadband Emissivity of the Sahara Desert Using ASTER
and MODIS Data," (in English), Earth Interactions, vol. 8, pp. -, 2004.
[60] A. M. Baldridge, S. J. Hook, C. I. Grove, and G. Rivera, "The ASTER Spectral Library
Version 2.0," Remote Sensing of Environment, vol. 114, no. 4, pp. 711-715, 2009.
[61] A. N. French, T. J. Schmugge, J. C. Ritchie, A. Hsu, F. Jacob, and K. Ogawa, "Detecting
land cover change at the Jornada Experimental Range, New Mexico with ASTER
emissivities," (in English), Remote Sensing of Environment, vol. 112, no. 4, pp. 1730-
1748, Apr 15 2008.
[62] M. Mira, E. Valor, R. Boluda, V. Caselles, and C. Coll, "Influence of soil water content
on the thermal infrared emissivity of bare soils: Implication for land surface temperature
determination," (in English), Journal of Geophysical Research-Earth Surface, vol. 112,
no. F4, pp. -, Oct 23 2007.
[63] A. N. French and A. Inamdar, "Land cover characterization for hydrological modelling
using thermal infrared emissivities," (in English), International Journal of Remote
Sensing, vol. 31, no. 14, pp. 3867-3883, 2010.
[64] A. A. Vandegriend and M. Owe, "On the Relationship between Thermal Emissivity and
the Normalized Difference Vegetation Index for Natural Surfaces," (in English),
International Journal of Remote Sensing, vol. 14, no. 6, pp. 1119-1131, Apr 1993.
[65] G. C. Hulley, S. J. Hook, and A. M. Baldridge, "ASTER land surface emissivity database
of California and Nevada," (in English), Geophysical Research Letters, vol. 35, no. 13,
pp. L13401, doi: 10.1029/2008gl034507, Jul 15 2008.
[66] G. C. Hulley and S. J. Hook, "Generating Consistent Land Surface Temperature and
Emissivity Products Between ASTER and MODIS Data for Earth Science Research," (in
English), Ieee Transactions on Geoscience and Remote Sensing, vol. 49, no. 4, pp. 1304-
1315, Apr 2011.
[67] U. NASS, "Field Crops: Usual Planting and Harvesting Dates," USDA National
Agricultural Statistics Service, Agriculural Handbook, no. 628, 2010.

37
[68] K. G. Laraby, J. R. Schott, and N. Raqueno, "Developing a confidence metric for the
Landsat land surface temperature product," in SPIE Defense+ Security, 2016, pp.
98400C-98400C-14: International Society for Optics and Photonics.
[69] M. Cook, J. R. Schott, J. Mandel, and N. Raqueno, "Development of an operational
calibration methodology for the landsat thermal data archive and initial testing of the
atmospheric compensation component of a land surface temperature (LST) product from
the archive," Remote Sensing, vol. 6, no. 11, pp. 11244-11266, 2014.
[70] J. A. Galve, C. Coll, V. Caselles, and E. Valor, "An atmospheric radiosounding database
for generating land surface temperature algorithms," (in English), Ieee Transactions on
Geoscience and Remote Sensing, vol. 46, no. 5, pp. 1547-1557, May 2008.
[71] R. Alley and M. Jentoft-Nilsen, "Algorithm Theoretical Basis Document for: Brightness
Temperature," Jet Propulsion Lab., California Inst. of Tech.; Pasadena, CA, United
States20060040979, 1999 1999.
[72] J. Kjaersgaard and R. Allen, "Loading the Landsat thermal band preprocessed using the
EROS LPGS system," Kimberly R&E Center, University of Idaho. Report, 2009.
[73] J. A. Augustine, J. J. DeLuisi, and C. N. Long, "SURFRAD - A national surface radiation
budget network for atmospheric research," (in English), Bulletin of the American
Meteorological Society, vol. 81, no. 10, pp. 2341-2357, 2000/10//undefined 2000.
[74] J. A. Augustine, G. B. Hodges, C. R. Cornwall, J. J. Michalsky, and C. I. Medina, "An
update on SURFRAD - The GCOS Surface Radiation budget network for the continental
United States," (in English), Journal of Atmospheric and Oceanic Technology, vol. 22,
no. 10, pp. 1460-1472, Oct 2005.
[75] K. C. Wang and S. L. Liang, "Evaluation of ASTER and MODIS land surface
temperature and emissivity products using long-term surface longwave radiation
observations at SURFRAD sites," (in English), Remote Sensing of Environment, vol. 113,
no. 7, pp. 1556-1565, Jul 15 2009.
[76] P. C. Guillevic et al., "Validation of Land Surface Temperature products derived from the
Visible Infrared Imaging Radiometer Suite (VIIRS) using ground-based and heritage
satellite measurements," (in English), Remote Sensing of Environment, vol. 154, pp. 19-
37, Nov 2014.
[77] J. A. Augustine and E. G. Dutton, "Variability of the surface radiation budget over the
United States from 1996 through 2011 from high-quality measurements," (in English),
Journal of Geophysical Research-Atmospheres, vol. 118, no. 1, pp. 43-53, 2013/01/16/
2013.
[78] P. C. Guillevic et al., "Validation of Land Surface Temperature products derived from the
Visible Infrared Imaging Radiometer Suite (VIIRS) using ground-based and heritage
satellite measurements," Remote Sensing of Environment, vol. 154, pp. 19-37,
2014/11/01/ 2014.
[79] P. Krishnan et al., "Comparison of in-situ, aircraft, and satellite land surface temperature
measurements over a NOAA Climate Reference Network site," Remote Sensing of
Environment, vol. 165, pp. 249-264, 2015.
[80] K. Ogawa, T. Schmugge, and S. Rokugawa, "Estimating broadband emissivity of arid
regions and its seasonal variations using thermal infrared remote sensing," (in English),
Ieee Transactions on Geoscience and Remote Sensing, vol. 46, no. 2, pp. 334-343, Feb
2008.

38
[81] Z. Zhu and C. E. Woodcock, "Object-based cloud and cloud shadow detection in Landsat
imagery," Remote Sensing of Environment, vol. 118, pp. 83-94, 2012/03/15/ 2012.
[82] H. Tonooka and F. D. Palluconi, "Validation of ASTER/TIR standard atmospheric
correction using water surfaces," (in English), Ieee Transactions on Geoscience and
Remote Sensing, vol. 43, no. 12, pp. 2769-2777, Dec 2005.
[83] Z. M. Wan, Y. L. Zhang, Q. C. Zhang, and Z. L. Li, "Validation of the land-surface
temperature products retrieved from Terra Moderate Resolution Imaging
Spectroradiometer data," (in English), Remote Sensing of Environment, vol. 83, no. 1-2,
pp. 163-180, Nov 2002.
[84] S. J. Hook et al., "Retrieval of lake bulk and skin temperatures using Along-Track
Scanning Radiometer (ATSR-2) data: A case study using Lake Tahoe, California," (in
English), Journal of Atmospheric and Oceanic Technology, vol. 20, no. 4, pp. 534-548,
Apr 2003.
[85] S. J. Hook, R. G. Vaughan, H. Tonooka, and S. G. Schladow, "Absolute radiometric in-
flight validation of mid infrared and thermal infrared data from ASTER and MODIS on
the Terra spacecraft using the Lake Tahoe, CA/NV, USA, automated validation site,"
IEEE Transactions on Geoscience and Remote Sensing, vol. 45, no. 6, p. 1798, 2007.
[86] G. C. Hulley, S. J. Hook, and P. Schneider, "Optimized split-window coefficients for
deriving surface temperatures from inland water bodies," (in English), Remote Sensing of
Environment, vol. 115, no. 12, pp. 3758-3769, Dec 15 2011.
[87] G. Hulley, N. Malakar, T. Hughes, T. Islam, and S. Hook, "MODIS MOD21 Land
Surface Temperature and Emissivity Algorithm Theoretical Basis Document," Jet
Propulsion Laboratory, California Institute of Technology, 2012.
[88] G. C. Hulley and S. J. Hook, "A radiance-based method for estimating uncertainties in
the Atmospheric Infrared Sounder (AIRS) land surface temperature product," (in
English), Journal of Geophysical Research-Atmospheres, vol. 117, Oct 26 2012.
[89] Z. Wan and Z. L. Li, "Radiance-based validation of the V5 MODIS land-surface
temperature product," (in English), International Journal of Remote Sensing, vol. 29, no.
17-18, pp. 5373-5395, 2008.
[90] S. M. Li, Y. Y. Yu, D. L. Sun, D. Tarpley, X. W. Zhan, and L. Chiu, "Evaluation of 10
year AQUA/MODIS land surface temperature with SURFRAD observations," (in
English), International Journal of Remote Sensing, vol. 35, no. 3, pp. 830-856, Feb 1
2014.
[91] P. Guillevic et al., "Land Surface Temperature product validation using NOAA's surface
climate observation networks—Scaling methodology for the Visible Infrared Imager
Radiometer Suite (VIIRS)," Remote Sensing of Environment, vol. 124, pp. 282-298,
2012.
[92] G. Hulley, S. Hook, W. Johnson, P. Guillevic, and N. Malakar, "Hyperspectral Thermal
Emission Spectrometer (HyTES) Level-2 Land Surface Temperature and Emissivity
Algorithm Theoretical Basis Document," Pasadena, CA: Jet Propulsion Laboratory,
National Aeronautics and Space Administration2016.
[93] M. Martin and F.-M. Göttsche, "Satellite LST Validation Report," KITGlobTemp_DEL-
12_i2r1, 12-Aug-2016 2016, issue Rev. 2.1.

39
Table 1: Instrument characteristics of the Landsat suite of thermal infrared sensors.

Satellite Sensor Period of time Thermal bands Resolution

Landsat 4 Thematic Mapper 07/1982 to Band 6:10.4 - 12.5 µm 120m

(TM) 12/1993
Landsat 5 Thematic Mapper 03/1984 to Band 6:10.4 - 12.5 µm 120m
(TM) 01/2013
Enhanced Band 6: 10.4 - 12.5
Landsat 7 Thematic Mapper Since 4/15/1999 60m
µm
Plus (ETM+)
Band 10: 10.6-
11.2µm
Landsat 8 Thermal Infrared 100m
Since 2/11/2013
Sensor (TIRS) Band 11: 11.5-
12.5µm

Table 2. The regression coefficients for spectral adjustment of ASTER GEDv3 for Landsat

sensors.

Sensor 𝑐𝑐13 𝑐𝑐14 𝑐𝑐

Landsat 4 0.3222 0.6498 0.0272
Landsat 5 -0.0723 1.0521 0.0195
Landsat 7 0.2147 0.7789 0.0058
Landsat 8 0.6820 0.2578 0.0584

Table 3. Emissivity validation of Landsat 4, 5, 7 and 8 at the Algodones and Kelso sand dune
sites. The Lab emissivities convolved with respective Landsat response function, and the Landsat
derived emissivity from the spectrally adjusted ASTER GEDv3 product show good agreement
with an average absolute error of 0.54%. The absolute error (|Lab Convolved-ASTER Adjusted|
× 100) is indicated by the symbol |δε|% in percentage.
Emissivity
Landsat
Sites Lab |δε| %
(ASTER adjusted)
Landsat 4 0.944 0.952 0.81
Landsat 5 0.956 0.956 0.02
Algodones
Landsat 7 0.947 0.953 0.59
Landsat 8 0.941 0.950 0.83
Kelso Landsat 4 0.944 0.951 0.70
40
 Landsat 5 0.954 0.956 0.22
Landsat 7 0.943 0.952 0.86
Landsat 8 0.937 0.947 1.03
Average 0.54

Table 4. SURFRAD sites used for LST validation over land with corresponding path and row for

Landsat scenes [https://www.esrl.noaa.gov/gmd/grad/surfrad/sitepage.html].

Code Name Latitude Longitude Land cover Path/Row

Bondville,
BND 40.052082°N -88.373032°W Cropland/Vegetation 23/32
Illinois
Goodwin Creek,
GWN 34.255087°N -89.873707°W Grassland 23/36
Mississippi
Penn. State
PSU Univ., 40.720287°N -77.930990°W Cropland/Vegetation 16/32
Pennsylvania
Sioux Falls,
SXF 43.734311°N -96.623365°W Grassland 29/30
South Dakota
Lake Tahoe, -120. 051577
TAH 39.134523 °N Water 43/33
CA/NV °W
SAL Salton Sea, CA 33.1630°N - 115.7864°W Water 39/37

Table 5. LST validation summary showing bias and RMSE in retrieved LST (K) for Landsat 5

(L5), Landsat 7 (L7) at selected SURFRAD sites. The number of scenes are indicated by N.

Sites Bondville, IL Goodwin Creek, MS Penn State, PA Sioux Falls, SD

Sensors Bias RMSE N Bias RMSE N Bias RMSE N Bias RMSE N
L5 1.3 2.5 80 0.5 2.0 109 0.2 2.3 59 0.7 2.1 60
L7 1.1 2.5 69 1.1 2.3 88 0.3 2.1 63 1.1 2.6 48
Both 10
1.2 2.5 149 0.8 2.2 197 0.3 2.2 122 0.9 2.3 8

41
Table 6. Landsat 5 and 7 LST validation summary for two inland water bodies in California/Nevada for

years 2000-2011. The number of scenes are indicated by N.

L5 L7 Both
Lake Tahoe, CA/NV Bias (K) -0.2 0.4 0.2
(39°N,120°W) RMSE (K) 0.6 0.8 0.7
N 57 122 179

L5 L7 Both
Salton Sea, CA Bias (K) -0.4 0.3 0.0
(33.3°N, 115.8°W)
RMSE (K) 0.7 0.5 0.6
N 36 48 84

Figure 1: Spectral response functions for ASTER bands 13 and 14, and thermal band 6 for

Landsat 4 (L4), Landsat 5 (L5), and Landsat 7 (L7), and band 10 and band 11 of Landsat 8 are

shown as L8a and L8b respectively.

42
Figure 2. Landsat 7 emissivity (TIR 10 µm-12.6 µm) cutout over the southwest USA derived

from a spectral adjustment and collocation of the ASTER GEDv3 product (see text for detail).

Emissivities vary from areas of low emissivity over sand dunes consisting of feldspars (A) and

quartz (E), high emissivity mafic rocks in the Mojave Desert (B), water (C), and vegetation (D).

43
1
2 Figure 3. Landsat 4, 5, 7 and 8 emissivity derived from ASTER GEDv3 bands 13 and 14 are

3 indicated by red symbols. The lab measured spectra at full resolution (black line), the lab spectra

4 convolved with the Landsat TIR band 6 (L4, L5 and L7) and band 10 (L8) are also indicated by

5 blue symbols. The higher Landsat derived ASTER emissivities at Kelso are most likely due to

6 the presence of sparse vegetation (desert grass) which raises the emissivity as reported by Hulley

7 et al. (2009).

44
 8

9 Figure 4. Landsat 7 emissivities (TIR 10 µm-12.6 µm) estimated from a spectral and vegetation

10 adjustment of ASTER GEDv3 emissivities for April and August 2002 over corn croplands near

45
11 Copeland, Kansas, USA (see details in text). Subplots (a) and b) show Landsat NDVI, c) and d)

12 the corresponding adjusted emissivities, e) Percentage emissivity change from the baseline

13 ASTER GED v3, and f) scene pixel counts for (e).

14

15

16 Figure 5: Landsat 5 LST uncertainty distribution plotted versus Total Column Water (TCW) and

17 simulated LST using the ASTER GEDv3 (left) and an emissivity classification (right) for a

18 selection of different rocks, soils and vegetation from the ASTER spectral library.

46
19

20
21 Figure 6. An example of the Land Surface Temperature (LST) retrieved from the Landsat 5 band
22 6 on 2 February, 2005 over the southwestern USA. The scene consists of a wide variety of land
23 cover classes such as water (Salton Sea at the center), sand dunes (Algodones Dunes), and
24 natural vegetation and managed agricultural systems in the Imperial Valley south of Salton Sea.
25

26
27
28
29

47
30
31 Figure 7. Location of LST Validation sites used in the study consisting of two inland water
32 bodies (a, and d), and four SURFRAD sites (b, c, e, and f) (Please refer to Table 4 for details).
33
34
35

48
36
37 Figure 8. Landsat 5 and 7 retrieved Land Surface Temperature (LST) matchups with ground-
38 based LST at four selected SURFRAD sites a) Bondville, IL (BND) b) Goodwin Creek, MS
39 (GWN) c) Penn State, PA (PSU) and d) Sioux Falls, SD (SXF). The data spans years 2000 to
40 2012, except for SXF (2003-2012). The numbers of observations are indicated in the legends.

49
41
42 Figure 9. Climatology of LST bias (Landsat-SURFRAD) shown with conventional box plot
43 diagrams at four SURFRAD sites a) Bondville, IL (BND) b) Goodwin Creek, MS (GWN) c)
44 Penn State, PA (PSU) and d) Sioux Falls, SD (SXF).The crop sites (a, and c) show two bias
45 peaks per year, which are likely related to crop growth, harvesting and regrowth of vegetation
46 while grasslands (a, and d) show higher biases during summer. The data spans years 2000 to
47 2012, except for SXF (2003-2012).

50
48
49 Figure 10. Distribution of LST accuracy for Landsat 5 and 7 at four SURFRAD sites indicated in
50 Figure 9.

51
52 Figure 11. Matchups of in situ derived lake surface temperature and Landsat 5 and 7 retrieved

53 LST over (a) Lake Tahoe, CA/NV and (b) Salton Sea, CA, from 2000-2011.

54

51
55
56 Figure 12. Distribution of Landsat LST accuracy for Landsat 5 and 7, at the two lake sites, Lake

57 Tahoe, CA/NV, and Salton Sea, CA.

52
58
59 Figure 13. a) Google Earth Image showing the footprints of L5/L7 TIR sensors compared with
60 the ground sensor’s field of view (circle) at the DRA site. b) HyTES airborne LST retrieval at ~6
61 m resolution at DRA shows large spatial heterogeneity in temperature, which can be as large as
62 ~15 K within different sensors field of view due to darker and hotter exposed gravel as shown in
63 the bottom panels (Indicated by labels c, and d, in Figure a). Images of bare exposed gravel/rock
64 surrounding desert rock, around the lightly colored area, d) image of darker and hotter gravel
65 area surrounding the site at location.

66

53

View publication stats