Analyzing Land Surface Temperature (LST) with Landsat 8 Data in Google Earth Engine

Note: this code is adapted from Muhammad Ridho’s excellent Medium tutorial. I have converted his JavaScript code to Python using the geemap package, with only minor changes. I am quoting his writeup almost verbatim.

This notebook is exclusively for educational purposes. All credit for the work goes to Muhammad.

Land Surface Temperature (LST) is a crucial environmental parameter with applications ranging from climate change monitoring to urban planning and agriculture management. This summary provides a concise overview of how to calculate and analyze LST using Landsat 8 imagery and the powerful toolset of Google Earth Engine (GEE).

Landsat 8, a satellite operated by the United States Geological Survey (USGS) and NASA, provides a wealth of multispectral data. Its thermal infrared band, known as the Thermal Infrared Sensor (TIRS), is particularly valuable for estimating LST. Google Earth Engine is a cloud-based platform for geospatial analysis, making it an ideal choice for processing Landsat data and deriving LST information.

Step 1: Accessing Google Earth Engine

1.1 Sign up for an Earth Engine account if you don’t have one already.

1.2 Load the ee and geemap libraries in Google Colab.

import ee
import geemap

1.3 Authenticate your Earth Engine account in Google Colab.

ee.Authenticate()

1.4 Initialize your session.

ee.Initialize(project='remotesensing-musa6500')

Step 2: Define the Area of Interest (AOI)

In this part of the script, an Area of Interest (AOI) is defined using the ee.FeatureCollection function. The AOI is obtained from the “FAO/GAUL_SIMPLIFIED_500m/2015/level1” dataset, specifically filtering for the administrative division (ADM1_NAME) named ‘Daerah Istimewa Yogyakarta.’ This effectively creates a region of interest for subsequent analysis.

aoi = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level1") \
        .filter(ee.Filter.eq('ADM1_NAME', 'Daerah Istimewa Yogyakarta'))

Step 3: Import Landsat 8 Data

3.1 Define the start and end dates.

These two variables, startDate and endDate, can be used in various geospatial analyses, particularly in Google Earth Engine, to specify a specific time range for data retrieval, processing, or visualization. In this case, the time period defined by these dates spans from June 1, 2023, to September 21, 2023.

from datetime import datetime

startDate = datetime(2023, 6, 1)
endDate = datetime(2023, 9, 21)

# Format dates into strings that Earth Engine expects ("YYYY-MM-DD")
startDate= startDate.strftime('%Y-%m-%d')
endDate_str = endDate.strftime('%Y-%m-%d')

3.2 Apply scaling factors

Applying a scale factor involves adjusting pixel values in satellite imagery to ensure they represent physical properties accurately. For example, it’s used to convert digital numbers into physical units, like radiance or reflectance, making the data more meaningful for analysis.

# Function to apply scaling factors.
def apply_scale_factors(image):
    # Scale and offset values for optical bands
    optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)

    # Scale and offset values for thermal bands
    thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0)

    # Add scaled bands to the original image
    return image.addBands(optical_bands, None, True) \
                .addBands(thermal_bands, None, True)

3.3 Apply cloud masking

Cloud masking is a process to identify and remove cloud-affected pixels from satellite imagery. It’s essential for ensuring that the analysis focuses on clear and cloud-free areas, minimizing the impact of clouds and their shadows on the results.

# Function to mask clouds and cloud shadows in Landsat 8 imagery
def cloud_mask(image):
    # Define cloud shadow and cloud bitmask (Bits 3 and 5)
    cloud_shadow_bit_mask = 1 << 3
    cloud_bit_mask = 1 << 5

    # Select the Quality Assessment (QA) band for pixel quality information
    qa = image.select('QA_PIXEL')

    # Create a binary mask to identify clear conditions (both cloud and cloud shadow bits set to 0)
    mask = qa.bitwiseAnd(cloud_shadow_bit_mask).eq(0) \
                .And(qa.bitwiseAnd(cloud_bit_mask).eq(0))

    # Update the original image, masking out cloud and cloud shadow-affected pixels
    return image.updateMask(mask)

3.4 Import Landsat 8 imagery

Importing images involves fetching satellite or remote sensing data into a geospatial analysis platform like Google Earth Engine. After data import, creating visualization involves defining how the data should be visually represented, and specifying parameters like color bands and brightness levels for effective interpretation and analysis.

# Import and preprocess Landsat 8 imagery
image = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2") \
            .filterBounds(aoi) \
            .filterDate(startDate, endDate) \
            .map(apply_scale_factors) \
            .map(cloud_mask) \
            .median() \
            .clip(aoi)

Step 4: Calculate NDVI, Proportion of Vegetation, and Emissivity

Step 4.1 Calculate NDVI

NDVI (Normalized Difference Vegetation Index) is a widely used vegetation index in remote sensing and environmental science. It quantifies the presence and health of vegetation based on the reflectance of visible and near-infrared light. NDVI values typically range from -1 to 1, with higher values indicating healthier or denser vegetation.

NDVI is calculated using the following formula:

NDVI=(NIR+Red)/(NIR−Red)​

  • NIR (Near-Infrared) is the reflectance in the near-infrared part of the electromagnetic spectrum (typically Landsat band 5 or similar).
  • Red is the reflectance in the red part of the spectrum (typically Landsat band 4 or similar).
# Calculate Normalized Difference Vegetation Index (NDVI)
ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI')

Step 4.2. Calculate Minimum and Maximum NDVI Values

The minimum (min) and maximum (max) values of the Normalized Difference Vegetation Index (NDVI) are used to calculate the proportion of vegetation cover (PV) in the region of interest.

Minimum NDVI Calculation (ndviMin):

  • The script uses the reduceRegion method to calculate the minimum NDVI value within the specified AOI.
  • The reducer parameter is set to ee.Reducer.min(), indicating that we want to find the minimum value.
  • The geometry parameter is set to aoi, which defines the area of interest.
  • The scale parameter sets the pixel size for the analysis to 30 meters.
  • The maxPixels parameter specifies the maximum number of pixels to include in the calculation (1e9 represents one billion pixels, a large number to cover extensive areas).

Maximum NDVI Calculation (ndviMax):

  • Similarly, this part of the script calculates the maximum NDVI value within the same AOI.
  • It uses ee.Reducer.max() as the reducer to find the maximum value.
# Calculate the minimum and maximum NDVI values within the AOI
ndvi_min = ndvi.reduceRegion(
    reducer=ee.Reducer.min(),
    geometry=aoi,
    scale=30,
    maxPixels=1e9
).get('NDVI').getInfo()

ndvi_max = ndvi.reduceRegion(
    reducer=ee.Reducer.max(),
    geometry=aoi,
    scale=30,
    maxPixels=1e9
).get('NDVI').getInfo()

Step 4.3. Calculate Proportion of Vegetation (PV) and Emissivity (EM)

The proportion of vegetation cover is calculated using the following formula:

Proportion of Vegetation (PV) = (NDVI — NDVI_min) / (NDVI_max — NDVI_min)

Where:

  • NDVI represents the pixel’s NDVI value.
  • NDVI_min is the minimum NDVI value specified.
  • NDVI_max is the maximum NDVI value specified.

The Proportion of Vegetation (PV) is a metric used to quantify the relative abundance of vegetation within a specified area by analyzing Normalized Difference Vegetation Index (NDVI) values. It provides valuable insights into land cover and ecosystem health, with higher PV values indicating a greater presence of vegetation. On the other hand, Emissivity (EM) is a critical parameter for accurate Land Surface Temperature (LST) calculations. In the provided code, EM is computed as a function of PV, reflecting how efficiently a surface emits thermal radiation. EM values near 1.0 are typical for natural surfaces like soil and vegetation, while lower values are often associated with water bodies or urban areas. Both PV and EM play essential roles in remote sensing and environmental studies, contributing to a more comprehensive understanding of land characteristics and thermal behavior.

# Convert NDVI_min and NDVI_max to ee.Number for operations
ndvi_min = ee.Number(ndvi_min)
ndvi_max = ee.Number(ndvi_max)

# Fraction of Vegetation (FV) Calculation
# Formula: ((NDVI - NDVI_min) / (NDVI_max - NDVI_min))^2
fv = ndvi.subtract(ndvi_min).divide(ndvi_max.subtract(ndvi_min)).pow(2).rename('FV')

# Emissivity Calculation
# Formula: 0.004 * FV + 0.986
em = fv.multiply(0.004).add(0.986).rename('EM')

Step 5: Estimate Land Surface Temperature (LST)

Select Thermal Band: In this section, the script selects the thermal band (Band 10) from the input image and renames it as ‘thermal’. This band contains thermal infrared information.

# Select Thermal Band (Band 10) and Rename It
thermal = image.select('ST_B10').rename('thermal')

The script calculates Land Surface Temperature (LST) using Planck’s Law. The formula is provided as a mathematical expression. It involves using the thermal band (TB) and emissivity (em) values to compute LST.

# Land Surface Temperature (LST) Calculation
# Formula: (TB / (1 + (λ * (TB / 1.438)) * log(em))) - 273.15
lst = thermal.expression(
    '(TB / (1 + (0.00115 * (TB / 1.438)) * log(em))) - 273.15',
    {
        'TB': thermal.select('thermal'),  # Thermal band as the temperature in Kelvin
        'em': em  # Emissivity
    }
).rename('LST Yogyakarta 2023')

Step 6: Map Results with Custom Title and Legend

Now, we can combine our results with custom visualizations to view them on the map. (This can be done step-by-step throughout the script or at the end.)

First, we initialize a map using the geemap package. We add our AOI to it as a black line using an empty palette and center the map on our AOI.

Map = geemap.Map() # initialize a map
Map.addLayer(aoi, {}, 'AOI - Yogyakarta')
Map.centerObject(aoi, 10) # center the map on the aoi

Next, we add our true color imagery to the map.

# Define visualization parameters for True Color imagery (bands 4, 3, and 2)
visualization = {
  'bands': ['SR_B4', 'SR_B3', 'SR_B2'],
  'min': 0.0,
  'max': 0.15,
}

# Add the processed image to the map with the specified visualization
Map.addLayer(image, visualization, 'True Color 432')

We add the NDVI layer, too.

# Define NDVI Visualization Parameters
ndviPalette = {
 'min': -1,
 'max': 1,
 'palette': ['blue', 'white', 'green']
}

Map.addLayer(ndvi, ndviPalette, 'NDVI Yogyakarta')

Finally, we add the LST layer with a custom pallete and visualize the map with a custom title and legend.

# Visualization parameters for LST
lst_vis_params = {
    'min': 18.47,
    'max': 42.86,
    'palette': [
        '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
        '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
        '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
        'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
        'ff0000', 'de0101', 'c21301', 'a71001', '911003'
    ]
}

# Add the LST layer to the map with custom visualization parameters
Map.addLayer(lst, lst_vis_params, 'Land Surface Temperature 2023')

# Define the color palette for the legend (updated to include the missing color '210300')
palette = [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003', '210300'
]

# Convert hex colors to the format required by the add_legend method
legend_colors = ['#' + color for color in palette]

# Calculate the legend keys based on the number of colors
# Assuming minLST and maxLST represent the range of your data
minLST, maxLST = 15, 45
legend_keys = [str(round(minLST + i * (maxLST - minLST) / (len(palette) - 1), 2)) for i in range(len(palette))]

# Add the custom legend to the map
Map.add_legend(keys=legend_keys, colors=legend_colors, position='bottomright')

# Set the center of the map (Example coordinates for Yogyakarta)
Map.set_center(110.3668, -7.8032, 10)

# Display the map
Map