gee-dan-2024-01.js

gistfile1.txt

// PRVI LST Calc v24
// There are multiple ways to calculate LST
// 1. Cook paper method (directly available from files ("ST_B10" - see below)
// 2. NDVI method (need to add)
// 3. Old method w/atm. coefficients (see below)
// Can we use the Lu and Ld to get LST the old way?
// How to get emissivity and water content?
// 4. Split-window method (will work with LS8?)

// Some functions:
// 1. Apply scaling factors for Surface Reflectance prods.
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.*').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  var emisBand = image.select("ST_EMIS").multiply(0.0001);
  var LdBand = image.select("ST_DRAD").multiply(0.001);
  var LuBand = image.select("ST_URAD").multiply(0.001);
  thermalBands = thermalBands.subtract(273.15)
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true)
              .addBands(emisBand, null, true)
              .addBands(LuBand, null, true)
              .addBands(LdBand, null, true);
}
// 2. Apply scaling factors for raw products
function applyScaleFactorsRaw(image) {
  var Band10 = image.select("B10").multiply(0.0003342).add(0.1);
  var Band11 = image.select("B11").multiply(0.0003342).add(0.1); //Does this band even work for LS8?
  // Bands 4 and 5 for NDVI (get here or from Surface Refl. files?)
  var Band4 = image.select("B4").multiply(0.010195).add(-0.1);
  var Band5 = image.select("B5").multiply(0.0062387).add(-0.1);
  return image.addBands(Band10, null, true)
              .addBands(Band11, null, true)
              .addBands(Band4, null, true)
              .addBands(Band5, null, true);
}
// 3. Mask clouds (all products)
function maskCloud(image) {
  // Bits 2, 3, and 4 are cirrus, cloud, and cloud shadow, respectively.
  // var cloudCirrusBitMask = (1 << 2); 
  var cloudsBitMask = (1 << 3);
  var cloudShadowBitMask = (1 << 4);
  // Get the pixel QA band.
  var qa = image.select('QA_PIXEL');
  // the flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                  // .and(qa.bitwiseAnd(cloudCirrusBitMask).eq(0)) 
                  .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

// PRVI (can restrict to give just one tile or filter by path/row)
var PRVI = ee.Geometry.Polygon([[-66.3, 17.30], [-65.2, 17.30], [-65.2, 18.45], [-66.3, 18.45],[-66.3, 17.30]]) ; 
// Surface Reflectance dataset
var dataset = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filterDate('2021-01-01', '2023-12-31')
//  .filterBounds(PRVI)
// Path and row filters
    .filterMetadata('WRS_PATH',"equals",5)
    .filterMetadata('WRS_ROW',"equals",47)
  .map(applyScaleFactors)
;
// Raw scene dataset
var dataset_raw = ee.ImageCollection('LANDSAT/LC08/C02/T1')
  .filterDate('2021-01-01', '2023-12-31')
//  .filterBounds(PRVI)
// Path and row filters
    .filterMetadata('WRS_PATH',"equals",5)
    .filterMetadata('WRS_ROW',"equals",47)
  .map(applyScaleFactorsRaw)
;

// Mask collection (also create unmasked version to check mask)
// var nomask = dataset;
var masked = dataset.map(maskCloud);
var masked_raw = dataset_raw.map(maskCloud);
print(masked);
print(masked_raw);

//var masked_raw_B10 = masked_raw.select('B10');
//print(masked_raw_B10)

// For calculations, use image collection, then reduce
// TOA radiance function (apply function to collection)
var addB10_rad_surf = function(coll,coll_sr) {
var B10_rad_surf = coll.expression("(((masked_B10 - Lu)/(e*t))-((1-e)/e)*Ld)",
    {masked_B10: coll.select('B10'),
    Lu: coll_sr.select('ST_URAD'),
    Ld: coll_sr.select('ST_DRAD'),
//     Lu: 4.12,  // Upwelling radiance
//     Ld: 6.08, // Downwelling radiance
     t: 0.50, // Transparency
     e: 0.95 // Emissivity
    }).rename('B10_SURF');
     return coll.addBands(B10_rad_surf);
};

// apply function to  collection
var B10_rad_surf = (masked_raw,masked).map(addB10_rad_surf);
print(B10_rad_surf);
















// Once calcs. are complete, select a single image
// Define list of images from collection here
// You can then select a single image if needed
var list = masked.toList(100, 0) ;
var list_raw = masked_raw.toList(100, 0) ;
print(list)
print(list_raw)

// Select single test image from list using index
// "13" is a good test image for Ponce, PR
var single_image = ee.Image(list.get(13));
var single_image_raw = ee.Image(list_raw.get(13));
// var single_image_LST = single_image.select("ST_B10");
// var single_image_EMIS = single_image.select("ST_EMIS");
// print(single_image_EMIS)

// Plot test image (seems ok)
// Map.addLayer(single_image, {bands: ['SR_B4', 'SR_B3', 'SR_B2'], max: 0.4}, 'masked');
// Map.addLayer(single_image_raw, {bands: ['B4', 'B3', 'B2'], max: 15000}, 'masked');











var LST = dataset.select("ST_B10");
print(LST)
var max_LST = LST.max();
var EMIS = dataset.select("ST_EMIS");
var mean_EMIS = EMIS.mean();

// Calculate LST from atm. coeffs. (old method)
// START HERE on 1/26/24
// Apply the mask to the image and display the result.
//var masked = cloudy_scene.updateMask(mask);
//Map.addLayer(masked, {bands: ['B4', 'B3', 'B2'], max: 0.4}, 'masked');

// Perform calculations on image2 (unmasked)
// Extract coefficients from imput image
//print(image_raw.getInfo());
var rad_mult = (image_raw.get('RADIANCE_MULT_BAND_10'));
print(rad_mult)
var rad_add = (image_raw.get('RADIANCE_ADD_BAND_10'));
print(rad_add)
var k1 = (image_raw.get('K1_CONSTANT_BAND_10'));
print(k1)
var k2 = (image_raw.get('K2_CONSTANT_BAND_10'));
print(k2)
var dt = (image_raw.get('DATE_ACQUIRED'));
print(dt)
var st = (image_raw.get('SCENE_CENTER_TIME'));
print(st)








// === add map =======================================================

// New palettes 
var palettes = require('users/gena/packages:palettes');
var palette = palettes.kovesi.rainbow_bgyr_35_85_c73[7];

//var LSTParams = {min:35, max:55,'palette':'000000, FF0000'};
var LSTParams = {min:30, max:60,palette:palette};
var EMISParams = {min:0.8, max:1.0,palette:palette};

// thermal viz:
Map.addLayer(max_LST,LSTParams,"max LST (DegC)");

//Map.addLayer(single_image_EMIS,EMISParams,"EMIS single image");
Map.addLayer(mean_EMIS,EMISParams,"EMIS mean");
Map.setCenter(-66.12,18.4, 11);