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

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. Go to the Google Earth Engine platform (https://earthengine.google.com/).

1.2. Sign in using your Google account or create one if you don’t have it.

Step 2: Define your AOI

2.1. Define the AOI (Area of Interest)

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

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.

2.2. Add the AOI to the Map and Set the Map View:

Map.addLayer(aoi, {}, 'AOI - Yogyakarta');
Map.centerObject(aoi, 10);

Here, the script adds the defined AOI to the map using Map.addLayer. The {} empty object is passed as the second argument, which means no special styling is applied to the AOI layer. The third argument, 'AOI - Yogyakarta', is the name that will be displayed for this layer on the map.

Next, Map.centerObject is used to center the map view on the AOI. The aoi is specified as the object to be centered and 10 is the zoom level.

Step 3: Import Landsat 8 Data

3.1. Define Start and End Date

var startDate = ‘2023–06–01’;
var endDate = ‘2023–09–21’;

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.

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.

// Applies scaling factors.
function applyScaleFactors(image) {
 // Scale and offset values for optical bands
 var opticalBands = image.select(‘SR_B.’).multiply(0.0000275).add(-0.2);
 
 // Scale and offset values for thermal bands
 var thermalBands = image.select(‘ST_B.*’).multiply(0.00341802).add(149.0);
 
 // Add scaled bands to the original image
 return image.addBands(opticalBands, null, true)
 .addBands(thermalBands, null, true);
}

3.3. Apply CloudMask

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

function cloudMask(image) {
  // Define cloud shadow and cloud bitmasks (Bits 3 and 5)
  var cloudShadowBitmask = (1 << 3);
  var cloudBitmask = (1 << 5);

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

  // Create a binary mask to identify clear conditions (both cloud and cloud shadow bits set to 0)
  var mask = qa.bitwiseAnd(cloudShadowBitmask).eq(0)
                .and(qa.bitwiseAnd(cloudBitmask).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
var image = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
              .filterBounds(aoi)
              .filterDate(startDate, endDate)
              .map(applyScaleFactors)
              .map(cloudMask)
              .median()
              .clip(aoi);

// Define visualization parameters for True Color imagery (bands 4, 3, and 2)
var 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');

Step 4: Calculate NDVI, Proportion of Vegetation, 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)
var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');

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

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

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.

// Calculate the minimum NDVI value within the AOI
var ndviMin = ee.Number(ndvi.reduceRegion({
  reducer   : ee.Reducer.min(),
  geometry  : aoi,
  scale     : 30,
  maxPixels : 1e9
}).values().get(0));

// Calculate the maximum NDVI value within the AOI
var ndviMax = ee.Number(ndvi.reduceRegion({
  reducer   : ee.Reducer.max(),
  geometry  : aoi,
  scale     : 30,
  maxPixels : 1e9
}).values().get(0));

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.

You can also print the minimum and maximum NDVI values using this script:

// Print the Minimum and Maximum NDVI Values
print("Minimum NDVI:", ndviMin);
print("Maximum NDVI:", ndviMax);

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.

// Fraction of Vegetation (FV) Calculation
// Formula: ((NDVI - NDVI_min) / (NDVI_max - NDVI_min))^2

// Calculate the Fraction of Vegetation (FV) using the NDVI values within the specified range.
// NDVI_min represents the minimum NDVI value, and NDVI_max represents the maximum NDVI value
var fv = ((ndvi.subtract(ndviMin)).divide(ndviMax.subtract(ndviMin)))
          .pow(ee.Number(2))
          .rename('FV');
        
  
// Emissivity Calculation
// Formula: 0.004 * FV + 0.986

// Calculate Land Surface Emissivity (EM) using the Fraction of Vegetation (FV).
// The 0.004 coefficient represents the emissivity variation due to vegetation,
// and the 0.986 represents the base emissivity for other surfaces.

var em = fv.multiply(ee.Number(0.004)).add(ee.Number(0.986)).rename('EM');

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.

Step 5: Land Surface Temperature Estimation

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

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.

// Now, lets calculate the land surface temperature (LST)

// Formula: (TB / (1 + (λ * (TB / 1.438)) * ln(em))) - 273.15
var lst = thermal.expression(
  '(TB / (1 + (0.00115 * (TB / 1.438)) * log(em))) - 273.15', {
    'TB': thermal.select('thermal'), // Select the thermal band (TB)
    'em': em // Assign emissivity (em)
  }).rename('LST Yogyakarta 2023');

// Add the LST Layer to the Map with Custom Visualization Parameters
Map.addLayer(lst, {
  min: 18.47, // Minimum LST value
  max: 42.86, // Maximum LST value
  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'
  ]}, 'Land Surface Temperature 2023');

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.

You can also access the LST values obtained using the Inspector tool.

Step 6. Title and Legend for LST Visualization

The final step involves adding a title and a legend to enhance the visualization of LST.

// Create a Legend for Land Surface Temperature (LST) Visualization
var minLST = 15; // Minimum LST value
var maxLST = 45; // Maximum LST value

// Create a panel for the legend with styling
var legend = ui.Panel({
  style: {
    position: 'bottom-right',
    padding: '8px 15px',
    backgroundColor: 'white'
  }
});

// Create a title for the legend
var legendTitle = ui.Label({
  value: 'Land Surface Temperature (°C)',
  style: {
    fontWeight: 'bold',
    fontSize: '20px',
    margin: '0 0 4px 0',
    padding: '0'
  }
});

// Add the legend title to the legend panel
legend.add(legendTitle);

// Define a color palette for the legend
var 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'
];

// Calculate the step value for the legend
var step = (maxLST - minLST) / (palette.length - 1);

// Loop through the palette and create legend entries
for (var i = 0; i < palette.length; i++) {
  // Create a color box for each legend entry
  var colorBox = ui.Label({
    style: {
      backgroundColor: '#' + palette[i],
      padding: '8px',
      margin: '0 0 8px 0',
      width: '42px'
    }
  });

  // Create a label with LST values for each legend entry
  var legendLabel = ui.Label({
    value: (minLST + i * step).toFixed(2),
    style: { margin: '0 0 10px 6px' }
  });

  // Create a panel to arrange color box and label horizontally
  var legendPanel = ui.Panel({
    widgets: [colorBox, legendLabel],
    layout: ui.Panel.Layout.Flow('horizontal')
  });

  // Add the legend entry panel to the main legend panel
  legend.add(legendPanel);
}

// Add the legend to the Google Earth Engine map
Map.add(legend);

// Create a Map Title
var mapTitle = ui.Panel({
  style: {
    position: 'top-center',
    padding: '20px 20px'
  }
});
var mapTitle2 = ui.Label({
  value: 'Land Surface Temperature - 2023 - Yogyakarta',
  style: {
    fontWeight: 'bold',
    fontSize: '30px',
    margin: '8px 8px 8px 8px',
    padding: '0'
  }
});
mapTitle.add(mapTitle2);

// Add the map title to the Google Earth Engine map
Map.add(mapTitle);

Visualize the title and legend