Saturday, January 4, 2025

GEE - Measuring Land Surface Temprature from Landsat 8

 The Javascript code below will calculate Land Surface Temperature (LST) from landsat 8 imagery.



Adopted from this medium post: Analyzing Land Surface Temperature (LST) with Landsat 8 Data in Google Earth Engine - https://medium.com/@ridhomuh002/analyzing-land-surface-temperature-lst-with-landsat-8-data-in-google-earth-engine-f4dd7ca28e70
// Load AOI or ROI
var aoi = ee.FeatureCollection("users/forum2k9/Karu"); // Karu polygon - change if you need different 
Map.addLayer(aoi, {}, 'AOI - Karu');
Map.centerObject(aoi, 10);

// Define data retrieval Start and End Date
var startDate = '2023-03-01';
var endDate = '2023-10-31';


// Applies scaling factors. Recalculating from K to °C
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).subtract(273.15); // I'm used to the °C, In case of K delete the subtract part
 // Add scaled bands to the original image
 return image.addBands(opticalBands, null, true)
 .addBands(thermalBands, null, true);
}

// Function to Mask Clouds and Cloud Shadows in Landsat 8 Imagery
function cloudMask(image) {
  // Define cloud shadow and cloud bitmasks
  var cloudShadowBitmask = (1 << 4);
  var cloudBitmask = (1 << 3);
  var cirrus = (1 << 2);             // I add also mask for cirrus. Suprisely use of Bit 6 (Clear of any clouds) is not working - anybody know why? 
  // 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))
                .and(qa.bitwiseAnd(cirrus).eq(0));
  // Update the original image, masking out cloud and cloud shadow-affected pixels
  return image.updateMask(mask);
}

// 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.2,
};

// Add the processed image to the map with the specified visualization
Map.addLayer(image, visualization, 'True Color 432');
print(image);
// Add the LST Layer to the Map with Custom Visualization Parameters
Map.addLayer(image.select('ST_B10'), {
  min: 19, // Minimum LST value
  max: 43, // 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');




// // 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 Karu')


// Adopted from: https://github.com/panmoravec/GEE/blob/main/LST%20Landsat8


Export.image.toDrive({
  image: image.select('ST_B10'),
  description: 'LST_2023',
  scale: 30,
  region: aoi,
  fileFormat: 'GeoTIFF',
  maxPixels: 3100100100,  // set bigger in case of need
});


// Another version for.... Land surface temperature (LST) Monitoring using Landsat imagery using Google Earth Engine

1) https://code.earthengine.google.com/525ebe138c6d2198a0fd464f7cfc967a 

var LLST = require('users/sofiaermida/landsat_smw_lst:modules/Landsat_LST.js')
var Lan8 = LLST.collection('L8', '2021-1-01', '2021-12-31', geometry, true).select('LST').map(function(image){
  return image.subtract(273.15).clip(geometry).set('system:time_start',image.get('system:time_start'))
})
print("chart");
// print(Lan8)
print(
  ui.Chart.image.series({
    imageCollection: Lan8,
    region: geometry,
    reducer: ee.Reducer.mean(),
    scale: 30,
    xProperty: 'system:time_start'
  }));


Map.centerObject(geometry,10)

Map.addLayer(geometry)

Map.addLayer(Lan8.mean().clip(geometry), imageVisParam4 ,'Dhaka City LST')


====================================================================================


2) https://code.earthengine.google.com/3fdeff0b82b6fc98022a11facd30d269 


// Land surface temperature
var landsat8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
             .filterBounds(ROI)
             .filterDate('2020-01-01','2022-12-31')
             //.filterMetadata('CLOUD_COVER', 'less_than', 1)
             .median()
             .clip(ROI)
print(landsat8)
Map.addLayer(landsat8, imageVisParam3, 'Satelliteimage with clouds')
Map.centerObject(landsat8,8)


// CLOUD MASK
//Bits 3 and 5 are cloud shadow and cloud, respectively
var  maskL8sr = function(col){
  var cloudShadowBitMask = (1 << 3)
  var cloudsBitMask = (1 << 5)
  
  // Get the piixel QA band
  var qa  = col.select('pixel_qa')
  
   var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
            .and(qa.bitwiseAnd(cloudsBitMask).eq(0))
  return col.updateMask(mask)
}

// visualization params
var vizParams = {
bands: ['B5', 'B6', 'B4'],
min: 0,
max: 4000,
gamma: [1, 0.9, 1.1]
};

var vizParams2 = {
bands: ['B4', 'B3', 'B2'],
min: 0,
max: 3000,
gamma: 1.4,
};

//load the collection
{
  var col = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
     .map(maskL8sr)
     .filterBounds(ROI)
     .filterDate('2020-01-01','2022-12-31')
     .mean()
     .clip(ROI)
}
print(col)
Map.addLayer(col, imageVisParam2, 'Cloud_Free')

// NDVI CALCULATION

var NDVI = col.normalizedDifference(['B5','B4'])
Map.addLayer(NDVI, {min: -1, max: 0.7, palette: ['blue','white','yellow','red']}, 'Created NDVI')
Map.centerObject(NDVI, 8)

print(NDVI)

// selecting Band10 which is the surface temperature
// thermal band 10
var thermalband  = col.select('B10').multiply(0.1)
Map.addLayer(thermalband, {min:294.388, max:300.514, palette:['yellow','green','blue']}, 'Thermal band')

//finding minimum and maximum NDVI

var min = ee.Number(NDVI.reduceRegion({   //reduce region is used to calculate a statistcal function for each region
reducer: ee.Reducer.min(),
geometry: ROI,
scale: 30,
maxPixels: 1e9
}).values().get(0));
print(min, 'min');
  
var max = ee.Number(NDVI.reduceRegion({
reducer: ee.Reducer.max(),
geometry: ROI,
scale: 30,
maxPixels: 1e9
}).values().get(0));
print(max, 'max')
  
//fractional vegetation
var fv = (NDVI.subtract(min).divide(max.subtract(min))).pow(ee.Number(2))
print(fv,'fv');
Map.addLayer(fv, {},'Fractional Vegetation');

//Emissivity
// you will have to find the emissivity of natural and town
 var a = ee.Number(0.004)
 var b = ee.Number(0.986)
 var EM  = fv.multiply(a).add(b).rename('EMM')
 Map.addLayer(EM, imageVisParam, 'Emissivity')
 
 //FINDING LAND SURFACE TEMPERATURE
 //LST in celsius degree bring -273.15
 //NB: In kelvin don't bring -273.15
 var LST = thermalband.expression(
   '(Tb/(1 + (0.00115*(Tb/1.438))*log(Ep)))-273.15',{
     'Tb':thermalband.select('B10'),
     'Ep':EM.select('EMM')
     
   }).rename('Land surface Temperature');
   var viz =  {min: 22.298638087799763, max:29.22115281845691, 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'
 ]}
   Map.addLayer(LST, viz, 'LST')
   
   
//TIMESSERIES


//print(image_drycollection)

//print(ui.Chart.image.series(LST, ROI, ee.Reducer.mean(), 1000, 'system:time_start'))

// Create the title label.
var title = ui.Label('Shama LandSurface Temperature Map Year of 2022 ');
title.style().set('position', 'top-center');
Map.add(title);

var title = ui.Label('Prepared by Kangah Desmond');
title.style().set('position', 'bottom-center');
Map.add(title);

// set position of panel
var legend = ui.Panel({
style: {
position: 'bottom-left',
padding: '8px 15px'
}
});

// Create the title label.
//var title = ui.Label('Dhaka City LandSurface Temperature Map Year of 2020 ');
//title.style().set('position', 'top-center');
//Map.add(title);
 
// Create legend title
var legendTitle = ui.Label({
value: 'LandSurface Temperature (degC)',
style: {
fontWeight: 'bold',
fontSize: '18px',
margin: '0 0 4px 0',
padding: '0'
}
});
 
// Add the title to the panel
legend.add(legendTitle);
 
// create the legend image
var lon = ee.Image.pixelLonLat().select('latitude');
var gradient = lon.multiply((viz.max-viz.min)/100.0).add(viz.min);
var legendImage = gradient.visualize(viz);
 
// create text on top of legend
var panel = ui.Panel({
widgets: [
ui.Label(viz['max'])
],
});
 
legend.add(panel);
 
// create thumbnail from the image
var thumbnail = ui.Thumbnail({
image: legendImage,
params: {bbox:'0,0,10,100', dimensions:'10x200'},
style: {padding: '1px', position: 'bottom-center'}
});
 
// add the thumbnail to the legend
legend.add(thumbnail);
 
// create text on top of legend
var panel = ui.Panel({
widgets: [
ui.Label(viz['min'])
],
});
 
legend.add(panel);
 
Map.add(legend);

// Create the title label.
//var title = ui.Label('Prepared By: Mijanur Rahman');
//title.style().set('position', 'bottom-right');
//Map.add(title);



Thank you for reading.

Wednesday, January 1, 2025

How to load shapefile into Google Earth Engine (GEE)

 To load a shapefile into Google Earth Engine (GEE), you can use the Code Editor, the Asset Manager, or the command line interface (CLI):

In this guide, I will use the 'Asset Manager' tab and access the shapefile from the Code Editor. Let get started...

Step 1: Go to code.earthengine.google.com and sign in. Click the “Assets” tab at the top left of the screen, and then click on the red “New” button. Select the “Shape files” option from the dropdown menu.


Step 2: In the pop-up window, click the red “Select” button and navigate to where your shapefile is stored. Select all of the files that make up your shapefile (.shp, .shx, .prj, .dbf, .cpg). You can hold the control key to select multiple files at once. Then click “Open” to add the selected files to the upload list. Update the Asset Name in the Asset ID section if needed. Then, click the “Upload” button at the bottom.



Step 3: Navigate to the 'Task' tab, there you see the submitted task in action, you don't have to do anything. Just watch the task will turn blue when the upload is complete.


Step 4: Click on the name of the name of the newly uploaded shapefile under the “Assets” tab to open its details windows where you will all meta data related to it including its "Table ID". The "Table ID" is what we will copy for use in the javascript code editor.



That is it!