[GIS] Calculating LST from Landsat 8 in Google Earth Engine

google-earth-enginejavascriptlandsatlandsat 8

I'm trying to write a code for Land Surface Temperature (LST) from Landsat 8 images in google earth engine. I'm using the code of this question as a guide and I succeeded at intermediate steps but I couldn't solve the previous issues of the script nor make it work.

Here's the link to the code:
https://code.earthengine.google.com/2b2c4e1ece6ac2a201f8e854c34cd577

Any suggestions/ideas for getting the final LST product?
My final goal is to obtain land surface temperature (from landsat 8), so if you know another way to get there that would be helpful too.

Here's my code:

//cloud mask
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

//vis 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)
.filterDate('2018-01-01','2018-12-31')
.filterBounds(geometry);
}
print(col, 'coleccion');

//median
{
var image = col.median();
print(image, 'image');
Map.addLayer(image, vizParams2);
}

// NDVI:
{
var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
print(ndvi,'ndvi');
Map.addLayer(ndvi, ndviParams, 'ndvi');
}

//

//select thermal band 10(with brightness tempereature), no BT calculation 
 var thermal= image.select('B10').multiply(1000);
 Map.addLayer(thermal, imageVisParam, 'thermal');


// find the min and max of NDVI
{
var min = ee.Number(ndvi.reduceRegion({
   reducer: ee.Reducer.min(),
   geometry: geometry,
   scale: 30,
   maxPixels: 1e9
   }).values().get(0));
print(min, 'min');
var max = ee.Number(ndvi.reduceRegion({
    reducer: ee.Reducer.max(),
   geometry: geometry,
   scale: 30,
   maxPixels: 1e9
   }).values().get(0));
print(max, 'max')
}

//fractional vegetation
{
var fv = ndvi.subtract(min).divide(max.subtract(min)).rename('FV'); 
print(fv, 'fv');
Map.addLayer(fv);
}

/////////////


  //Emissivity

  var a= ee.Number(0.004);
  var b= ee.Number(0.986);
  var EM=fv.multiply(a).add(b).rename('EMM');
  Map.addLayer(EM, imageVisParam2,'EMM');


  //LST c,d,f, p1, p2, p3 are assigned variables to write equaton easily
  var c= ee.Number(1);
  var d= ee.Number(0.00115);
  var f= ee.Number(1.4388);


var p1= ee.Number(thermal.multiply(d).divide(f));
var p2= ee.Number(Math.log(EM));
var p3= ee.Number((p1.multiply(p2)).add(c));


  var LST= (thermal.divide(p3)).rename('LST');

  var LSTimage = ee.Image(LST)

  Map.addLayer(LSTimage, {min: 0, max: 350, palette: ['FF0000', 
  '00FF00']},'LST');

Best Answer

In the last answer, he forgot to do the image reduction. For a better result use this code:

//cloud mask
function maskL8sr(col) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = col.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return col.updateMask(mask);
}

//vis 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)
.filterDate('2018-01-01','2018-12-31')
.filterBounds(geometry);
}
print(col, 'coleccion');

//imagen reduction
{
var image = col.median();
print(image, 'image');
Map.addLayer(image, vizParams2);
}


//median
{
var ndvi = image.normalizedDifference(['B5', 
'B4']).rename('NDVI');
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 
'green']};
print(ndvi,'ndvi');
Map.addLayer(ndvi, ndviParams, 'ndvi');
}

//select thermal band 10(with brightness tempereature), no calculation 
var thermal= image.select('B10').multiply(0.1);
var b10Params = {min: 291.918, max: 302.382, palette: ['blue', 
'white', 'green']};
Map.addLayer(thermal, b10Params, 'thermal');

// find the min and max of NDVI
{
var min = ee.Number(ndvi.reduceRegion({
reducer: ee.Reducer.min(),
geometry: geometry,
scale: 30,
maxPixels: 1e9
}).values().get(0));
print(min, 'min');
var max = ee.Number(ndvi.reduceRegion({
reducer: ee.Reducer.max(),
geometry: geometry,
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)).rename('FV'); 
print(fv, 'fv');
Map.addLayer(fv);
}

//Emissivity

var a= ee.Number(0.004);
var b= ee.Number(0.986);
var EM=fv.multiply(a).add(b).rename('EMM');
var imageVisParam3 = {min: 0.9865619146722164, max:0.989699971371314};
Map.addLayer(EM, imageVisParam3,'EMM');

//LST in Celsius Degree bring -273.15
//NB: In Kelvin don't bring -273.15
var LST = thermal.expression(
'(Tb/(1 + (0.00115* (Tb / 1.438))*log(Ep)))-273.15', {
 'Tb': thermal.select('B10'),
'Ep': EM.select('EMM')
}).rename('LST');
Map.addLayer(LST, {min: 20.569706944223423, max:29.328077233404645, 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'
 ]},'LST');
Related Question