Google Earth Engine Code

[Link]
66946d98f83c378d6a5616bbeffdc899

// ==============================================================================

// STEP 1: Define Study Area (BMC Boundaries) and Load Satellite Data

// ==============================================================================

// 1. Load the Global Administrative Boundaries dataset for Maharashtra

var maharashtra = [Link]("projects/fourth-vehicle-452506-f2/assets/shape1");

// 2. Set this exact combined boundary as your new Region of Interest (ROI)

var roi = [Link]();

// 3. Center the map and draw a black outline

[Link](roi, 11);

[Link]([Link]().paint(roi, 0, 2), {palette: 'black'}, 'BMC Boundary');

function maskL8(image){

var qa = [Link]('QA_PIXEL');

var cloud = [Link](1 << 3).eq(0);

var cloudShadow = [Link](1 << 4).eq(0);

return [Link](cloud).updateMask(cloudShadow);

// 4. Fetch the raw collection FIRST so we can read its metadata (NOW UPDATED TO 2025)

var landsatCollection = [Link]('LANDSAT/LC08/C02/T1_L2')

.filterBounds(roi)

.filterDate('2025-01-01', '2025-05-31')
 .map(maskL8);

// Bumped to 15% to ensure a clean mosaic for 2025

// 5. Print all the individual images and their data to the Console

print('--- SATELLITE IMAGE METADATA ---');

print('Images used for this map:', landsatCollection);

// 6. Smash them together into a single, clear median image for the classification

// 6. Define the function to calculate NDVI, NDBI, and NDWI

var addIndices = function(image){

var ndvi = [Link](['SR_B5','SR_B4']).rename('NDVI');

var ndbi = [Link](['SR_B6','SR_B5']).rename('NDBI');

var ndwi = [Link](['SR_B3','SR_B5']).rename('NDWI');

return [Link]([ndvi, ndbi, ndwi]);

};

// 7. Smash them together, scale the data to true reflectance, and add indices!

var landsat = landsatCollection

.median()

.clip(roi)

.multiply(0.0000275).add(-0.2); // Official Landsat 8 scaling

landsat = addIndices(landsat);

// ==============================================================================

// EXTRA: Calculate the Total Area of the Study Boundary (ROI)

// ==============================================================================

// 1. Earth Engine calculates area in square meters by default.

// 2. We divide by 1,000,000 (1e6) to convert it to Square Kilometers.

var totalAreaSqKm = [Link]().divide(1e6);

print('--- TOTAL STUDY AREA ---');

print('Total BMC Boundary Area (Sq Km):', totalAreaSqKm);

// ==============================================================================

// STEP 2: Train the Classifier (70/30 Split) and Build the LULC Map

// ==============================================================================

// 1. Merge all your manual pins together

var allPoints = [Link](forest)

.merge(water)

.merge(vegetation)

.merge(agriculture);

// 2. Add a random number to each pin so we can shuffle them fairly

var withRandom = [Link]('random_number', 42);

// 3. Split the data! 70% goes to Training, 30% is hidden for Validation

var trainingSample = [Link]([Link]('random_number', 0.7));

var validationSample = [Link]([Link]('random_number', 0.7));

var bands = ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7', 'NDVI', 'NDBI', 'NDWI'];

// 4. Extract the satellite data specifically for the 70% TRAINING set

var trainingData = [Link](bands).sampleRegions({

collection: trainingSample,

properties: ['class'],

scale: 30

});

// 5. TRAIN THE MODEL! (Using 100 decision trees for better accuracy)

var classifier = [Link]({

numberOfTrees: 300,
 variablesPerSplit: 4,

minLeafPopulation: 1,

bagFraction: 0.7

}).train({

features: trainingData,

classProperty: 'class',

inputProperties: bands

});

// 6. Generate your final map using the newly trained model

var myCustomMap = [Link](bands).classify(classifier).rename('classification');

// ==============================================================================

// STEP 2B: TEST THE MODEL (The 0.90+ Kappa Check on Hidden Data)

// ==============================================================================

// 1. Extract the satellite data for the hidden 30% VALIDATION set

var validationData = [Link](bands).sampleRegions({

collection: validationSample,

properties: ['class'],

scale: 30

});

// 2. Force the trained model to guess the classes of your hidden data

var validated = [Link](classifier);

// 3. Generate the Confusion Matrix to see how well it scored!

var validationMatrix = [Link]('class', 'classification');

print('--- MY CUSTOM MODEL ACCURACY (70/30 Split) ---');

print('Validation Matrix:', validationMatrix);

print('Overall Accuracy:', [Link]());

print('Kappa Coefficient:', [Link]());

// ADDED: Print the detailed class-by-class accuracy!

// Output will be a list in order of your classes: [Urban, Forest, Water, Veg, Agri]

print('Producer\'s Accuracy (Omission Error):', [Link]());

print('User\'s Accuracy (Commission Error):', [Link]());

// ==============================================================================

// STEP 2C: HISTORICAL LULC CLASSIFICATION (1990, 2000, 2015)

// ==============================================================================

print('Generating Historical Maps & Applying Reverse Sprawl Filter...');

// 1. Master Function to fetch Harmonized Landsat data

var getHarmonizedLandsat = function(year) {

var startDate = [Link](year, 1, 1);

var endDate = [Link](year, 5, 31);

var image;

var bandsToSelect;

var standardNames = ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'];

if (year < 2013) {

image = [Link]('LANDSAT/LT05/C02/T1_L2')

.filterBounds(roi).filterDate(startDate, endDate).median();

bandsToSelect = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'];

} else {

image = [Link]('LANDSAT/LC08/C02/T1_L2')

.filterBounds(roi).filterDate(startDate, endDate).median();

bandsToSelect = ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'];

}
 image = [Link](bandsToSelect, standardNames).clip(roi).multiply(0.0000275).add(-0.2);

var ndvi = [Link](['NIR', 'Red']).rename('NDVI');

var ndbi = [Link](['SWIR1', 'NIR']).rename('NDBI');

var ndwi = [Link](['Green', 'SWIR1']).rename('NDWI');

return [Link]([ndvi, ndbi, ndwi]);

};

// 2. Train and Classify

var generateHistoricalMap = function(year, trainingPins) {

var historicalLandsat = getHarmonizedLandsat(year);

var predictionBands = ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'NDVI', 'NDBI', 'NDWI'];

var trainingData = [Link](predictionBands).sampleRegions({

collection: trainingPins, properties: ['class'], scale: 30

});

var historicalClassifier = [Link]({

numberOfTrees: 100, variablesPerSplit: 3

}).train({

features: trainingData, classProperty: 'class', inputProperties: predictionBands

});

return [Link](predictionBands).classify(historicalClassifier).rename('classification');

};

// 3. GENERATE RAW MAPS

var raw1990 = generateHistoricalMap(1990, allPoints);

var raw2000 = generateHistoricalMap(2000, allPoints);

var raw2015 = generateHistoricalMap(2015, allPoints);

// 4. THE REVERSE SPRAWL FILTER (Ensures 1990 < 2000 < 2015 < 2016)

// First, fetch the strict 2016 Dynamic World Baseline

var dw16_baseline = [Link]("GOOGLE/DYNAMICWORLD/V1")

.filterBounds(roi).filterDate('2016-01-01', '2016-12-31').select('label').mode().clip(roi)

.remap([6, 1, 0, 2, 3, 5, 4, 7, 8], [1, 2, 3, 4, 4, 4, 5, 5, 5]).rename('class');

// Filter backwards: If it wasn't urban in the future, force it to match nature in the past

var map2015 = [Link]([Link](1).and(dw16_baseline.neq(1)), dw16_baseline);

var map2000 = [Link]([Link](1).and([Link](1)), map2015);

var map1990 = [Link]([Link](1).and([Link](1)), map2000);

var histPalette = ['red', 'green', 'blue', 'pink', 'yellow'];

[Link](map1990, {min: 1, max: 5, palette: histPalette}, 'Corrected LULC 1990', false);

[Link](map2000, {min: 1, max: 5, palette: histPalette}, 'Corrected LULC 2000', false);

[Link](map2015, {min: 1, max: 5, palette: histPalette}, 'Corrected LULC 2015', false);

// ==============================================================================

// STEP 2D: HISTORICAL LULC AREA STATISTICS (Horizontal Sq Km Layout)

// ==============================================================================

print('Calculating Filtered Historical Area Statistics (1990, 2000, 2015)...');

var getMapAreasSqKm = function(lulcImage) {

var areaImage = [Link]().divide(1e6).addBands([Link]('class'));

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'class'}),

geometry: roi, scale: 30, maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

var classAreas = [Link](function(item) {

 var d = [Link](item);

var classString = [Link]([Link]('class')).toInt().format('%d');

return [Link]([classString, [Link]('sum')]);

});

var defaultDict = [Link]({'1':0, '2':0, '3':0, '4':0, '5':0});

return [Link]([Link]([Link]()), true);

};

// Measure the newly filtered maps!

var area90 = getMapAreasSqKm(map1990);

var area00 = getMapAreasSqKm(map2000);

var area15 = getMapAreasSqKm(map2015);

var histClassNames = [Link](['1. Urban', '2. Forest', '3. Water', '4. Vegetation', '5. Agriculture']);

var histClassValues = [Link]([1, 2, 3, 4, 5]);

var table2Rows = [Link](function(c) {

var cNum = [Link](c);

var cKey = [Link]('%d');

var name = [Link]([Link](1));

var a90 = [Link]([Link](cKey));

var a00 = [Link]([Link](cKey));

var a15 = [Link]([Link](cKey));

var diff90_00 = [Link](a90);

var diff00_15 = [Link](a00);

var diff_Total_90_to_15 = [Link](a90);

return [Link](null, {
 'LULC_Class': name,

'Area_1990': [Link]('%.2f'),

'Area_2000': [Link]('%.2f'),

'Area_2015': [Link]('%.2f'),

'Change_90_to_00': diff90_00.format('%.2f').cat(' sq km'),

'Change_00_to_15': diff00_15.format('%.2f').cat(' sq km'),

'TOTAL_CHANGE_90_to_15': diff_Total_90_to_15.format('%.2f').cat(' sq km')

});

});

var paperTable2 = [Link]([Link](table2Rows), 'LULC_Class')

.setChartType('Table')

.setOptions({

title: 'Filtered Historical Land Cover Change in Sq Km (1990 - 2015)',

allowHtml: true, pageSize: 5

});

print('--- HISTORICAL LULC STATISTICS ---');

print(paperTable2);

// ==============================================================================

// EXTRA: Calculate and Graph the Area of Your Custom 2025 Map

// ==============================================================================

var areaImage = [Link]().divide(1e6);

var areaByClass = [Link](myCustomMap).reduceRegion({

reducer: [Link]().group({

groupField: 1,

groupName: 'class_number'

}),

geometry: roi,

scale: 30,
 maxPixels: 1e10

});

var classAreas = [Link]([Link]('groups')).map(function(item) {

var areaDict = [Link](item);

var classNumber = [Link]([Link]('class_number'));

var areaSqKm = [Link]([Link]('sum'));

var classNames = [Link](['Urban', 'Forest', 'Water', 'Vegetation', 'Agriculture']);

var name = [Link]([Link](1));

return [Link](null, {'Land Cover': name, 'Area (Sq Km)': areaSqKm});

});

var customAreaChart = [Link]([Link](classAreas), 'Land Cover',

'Area (Sq Km)')

.setChartType('ColumnChart')

.setOptions({

title: 'Land Cover Area Breakdown (My Custom 2025 Map)',

hAxis: {title: 'Land Cover Class'},

vAxis: {title: 'Total Area (Square Kilometers)'},

colors: ['#8A2BE2'], // Using a distinct color for your custom 2025 map

legend: {position: 'none'}

});

print('--- MY CUSTOM MAP AREA (2025) ---');

print(customAreaChart);

// ==============================================================================

// GRAPH: USER'S AND PRODUCER'S ACCURACY BY CLASS (SKIPPING CLASS 0)

// ==============================================================================
// 1. Flatten the complex matrices into simple 1D lists

var pAccList = [Link]([Link]().project([0]).toList());

var uAccList = [Link]([Link]().project([1]).toList());

// 2. Define your exact class names (5 classes)

var classNames = [Link](['Urban', 'Forest', 'Water', 'Vegetation', 'Agriculture']);

// 3. Loop through indices 1 to 5 (Completely ignoring 0!)

var indices = [Link](1, 5);

var accuracyFeatures = [Link](function(index) {

var n = [Link](index);

// To match index 1 to the 1st name in our list, we subtract 1

var name = [Link]([Link](1));

// Pull the exact percentage using the actual class numbers (1, 2, 3, 4, 5)

var pAcc = [Link]([Link](n)).multiply(100);

var uAcc = [Link]([Link](n)).multiply(100);

// Bundle them into a Feature

return [Link](null, {

'Class': name,

'Producer Accuracy (%)': pAcc,

'User Accuracy (%)': uAcc

});

});

// 4. Build the Professional Bar Chart

var accuracyChart = [Link]([Link](accuracyFeatures), 'Class',
['Producer Accuracy (%)', 'User Accuracy (%)'])

.setChartType('ColumnChart')

.setOptions({

title: 'Model Reliability vs. Completeness (By Class)',

hAxis: {title: 'Land Cover Class'},

vAxis: {

title: 'Accuracy Percentage (%)',

minValue: 0,

maxValue: 100

},

colors: ['#1f77b4', '#ff7f0e'],

legend: {position: 'bottom'}

});

print('--- ACCURACY GRAPH ---');

print(accuracyChart);

// ==============================================================================

// STEP 2E: HISTORICAL LAND SURFACE TEMPERATURE (1990, 2000, 2015)

// ==============================================================================

// ==============================================================================

// STEP 2E: HISTORICAL LAND SURFACE TEMPERATURE (Gap-Filled!)

// ==============================================================================

print('Calculating Historical LST Trends...');

// 1. Master Function to safely extract LST and fill missing data gaps

var getHistoricalLST = function(year) {

var startDate = [Link](year, 1, 1);

var endDate = [Link](year, 12, 31); // ← FULL YEAR to maximize coverage

var lSat;

var stBand;

if (year < 2013) {

// Merge ALL Landsat 5 tiers across the FULL YEAR

var t1_5 = [Link]('LANDSAT/LT05/C02/T1_L2');

var t2_5 = [Link]('LANDSAT/LT05/C02/T2_L2');

lSat = t1_5.merge(t2_5)

.filterBounds(roi)

.filterDate(startDate, endDate)

// Apply pixel-level QA cloud mask BEFORE median

.map(function(img) {

var qa = [Link]('QA_PIXEL');

var cloudFree = [Link](1 << 3).eq(0) // No clouds

.and([Link](1 << 4).eq(0)); // No cloud shadows

return [Link](cloudFree);

})

.median();

stBand = 'ST_B6';

} else {

var t1_8 = [Link]('LANDSAT/LC08/C02/T1_L2');

var t2_8 = [Link]('LANDSAT/LC08/C02/T2_L2');

lSat = t1_8.merge(t2_8)

.filterBounds(roi)

.filterDate(startDate, endDate)

.map(function(img) {

var qa = [Link]('QA_PIXEL');

var cloudFree = [Link](1 << 3).eq(0)

 .and([Link](1 << 4).eq(0));

return [Link](cloudFree);

})

.median();

stBand = 'ST_B10';

// Convert to Celsius

var lst = [Link](stBand)

.multiply(0.00341802).add(149.0).subtract(273.15);

// --- GAP FILL: Use focal_mean to interpolate missing pixels ---

// First pass: small radius to fill thin gaps

var filled1 = lst.focal_mean({radius: 300, units: 'meters', iterations: 3});

// Second pass: larger radius for bigger gaps (e.g. Borivali white patches)

var filled2 = filled1.focal_mean({radius: 1000, units: 'meters', iterations: 2});

// Use original where valid, filled where missing

var gapFilled = [Link](filled1).unmask(filled2);

// Mumbai realistic temp range: 18°C to 55°C (wide enough for all seasons)

return gapFilled

.updateMask([Link](18).and([Link](55)))

.rename('LST_' + year)

.clip(roi);

};

// 2. Extract the actual thermal maps for our three decades

var lst90 = getHistoricalLST(1990);

var lst00 = getHistoricalLST(2000);

var lst15 = getHistoricalLST(2015);

// 3. Add them to the Map Viewer (Turned off by default)

//var thermalPalette = ['blue', 'cyan', 'green', 'yellow', 'orange', 'red'];

//[Link](lst90, {min: 25, max: 45, palette: thermalPalette}, 'Historical LST 1990', false);

//[Link](lst00, {min: 25, max: 45, palette: thermalPalette}, 'Historical LST 2000', false);

//[Link](lst15, {min: 25, max: 45, palette: thermalPalette}, 'Historical LST 2015', false);

// 4. Helper Function: Calculate the Mean Temperature for a specific LULC class

// (This safely ties the temperature of that year to the specific LULC map of that year!)

var extractMeanLST = function(lstMap, lulcMap, classNum) {

var temp = [Link]([Link](classNum)).reduceRegion({

reducer: [Link](),

geometry: roi,

scale: 30,

maxPixels: 1e10

});

return [Link]().get(0);

};

// 5. Calculate the exact temperature for all 5 classes across the 3 historical years

var u90 = extractMeanLST(lst90, map1990, 1); var u00 = extractMeanLST(lst00, map2000, 1); var u15
= extractMeanLST(lst15, map2015, 1);

var f90 = extractMeanLST(lst90, map1990, 2); var f00 = extractMeanLST(lst00, map2000, 2); var f15 =
extractMeanLST(lst15, map2015, 2);

var w90 = extractMeanLST(lst90, map1990, 3); var w00 = extractMeanLST(lst00, map2000, 3); var
w15 = extractMeanLST(lst15, map2015, 3);

var v90 = extractMeanLST(lst90, map1990, 4); var v00 = extractMeanLST(lst00, map2000, 4); var v15
= extractMeanLST(lst15, map2015, 4);

var a90 = extractMeanLST(lst90, map1990, 5); var a00 = extractMeanLST(lst00, map2000, 5); var a15
= extractMeanLST(lst15, map2015, 5);
// 6. Build the Grouped Column Chart

var histTempChart = [Link]([Link]([

[Link](null, {'Class': '1. Urban', '1990 (°C)': u90, '2000 (°C)': u00, '2015 (°C)': u15}),

[Link](null, {'Class': '2. Forest', '1990 (°C)': f90, '2000 (°C)': f00, '2015 (°C)': f15}),

[Link](null, {'Class': '3. Water', '1990 (°C)': w90, '2000 (°C)': w00, '2015 (°C)': w15}),

[Link](null, {'Class': '4. Vegetation', '1990 (°C)': v90, '2000 (°C)': v00, '2015 (°C)': v15}),

[Link](null, {'Class': '5. Agriculture', '1990 (°C)': a90, '2000 (°C)': a00, '2015 (°C)': a15})

]), 'Class', ['1990 (°C)', '2000 (°C)', '2015 (°C)'])

.setChartType('ColumnChart')

.setOptions({

title: 'Historical Land Surface Temperature by Class (1990, 2000, 2015)',

hAxis: {title: 'Land Cover Class', textStyle: {bold: true}},

vAxis: {title: 'Mean LST (°C)', textStyle: {bold: true}},

// Visualizing the heat progression over time: Cool Blue (1990) -> Orange (2000) -> Hot Red (2015)

colors: ['#4682B4', '#FFA500', '#DC143C'],

legend: {position: 'top'}

});

print('--- HISTORICAL LST TREND ---');

print(histTempChart);

// ==============================================================================

// STEP 3: Dynamic World Time-Series & Multi-Year Accuracy Validation

// ==============================================================================

var getDynamicWorld = function(year) {

var y = [Link](year);

var startDate = [Link](y, 1, 1);

var endDate = [Link](y, 5, 31);

var dw = [Link]("GOOGLE/DYNAMICWORLD/V1")

.filterBounds(roi)
 .filterDate(startDate, endDate)

.select('label')

.mode()

.clip(roi);

return [Link](

[6, 1, 0, 2, 3, 5, 4, 7, 8],

[1, 2, 3, 4, 4, 4, 5, 5, 5]

).rename('dw_reference_class');

};

var dw2016 = getDynamicWorld(2016);

var dw2020 = getDynamicWorld(2020);

var dw2025 = getDynamicWorld(2025);

// --- GRADE GOOGLE DYNAMIC WORLD ---

// --- GRADE GOOGLE DYNAMIC WORLD ---

var groundTruth = allPoints;

var gradeDynamicWorld = function(dwImage, yearLabel) {

var dwTested = [Link]({

collection: groundTruth,

properties: ['class'],

scale: 10

});

var dwMatrix = [Link]('class', 'dw_reference_class');

print('--- GOOGLE DYNAMIC WORLD (' + yearLabel + ') ---');

print('Overall Accuracy:', [Link]());

 print('Kappa Coefficient:', [Link]());

// ADDED: Print the detailed class-by-class accuracy for Dynamic World!

print('Producer\'s Accuracy:', [Link]());

print('User\'s Accuracy:', [Link]());

};

gradeDynamicWorld(dw2016, '2016');

gradeDynamicWorld(dw2020, '2020');

gradeDynamicWorld(dw2025, '2025');

// ==============================================================================

// STEP 3B: Cross-Validation (My Custom Model vs. Dynamic World 2025)

// ==============================================================================

var combinedImage = [Link](dw2025); // Now checking against DW 2025!

var crossValPoints = [Link]({

region: roi,

points: 1000,

seed: 42

});

var crossValData = [Link]({

collection: crossValPoints,

scale: 30

});

var crossMatrix = [Link]({

actual: 'dw_reference_class',

predicted: 'classification'
});

print('--- CROSS-VALIDATION: MY MODEL vs. DYNAMIC WORLD (2025) ---');

print('Cross-Validation Matrix:', crossMatrix);

print('Cross-Validation Accuracy:', [Link]());

print('Cross-Validation Kappa:', [Link]());

// ==============================================================================

// GRAPH: OVERALL ACCURACY COMPARISON (Dynamic World vs. My Custom Model)

// ==============================================================================

// 1. Get the accuracy of your Custom Model (Calculated in Step 2B)

var myAcc = [Link]().multiply(100);

// 2. Create a quick helper function to extract the Overall Accuracy % for Dynamic World

var getDwAccuracy = function(dwImage) {

var dwTested = [Link]({

collection: allPoints, // Your hand-drawn ground truth pins

properties: ['class'],

scale: 10

});

return [Link]('class', 'dw_reference_class').accuracy().multiply(100);

};

// Extract the exact accuracy percentages for the historical maps

var acc16 = getDwAccuracy(dw2016);

var acc20 = getDwAccuracy(dw2020);

var acc25 = getDwAccuracy(dw2025);

// 3. Bundle them together into a dataset for the Chart

var accuracyComparison = [Link]([

 [Link](null, {'Map Version': 'DW 2016', 'Overall Accuracy (%)': acc16}),

[Link](null, {'Map Version': 'DW 2020', 'Overall Accuracy (%)': acc20}),

[Link](null, {'Map Version': 'DW 2025', 'Overall Accuracy (%)': acc25}),

[Link](null, {'Map Version': 'My Custom Map', 'Overall Accuracy (%)': myAcc})

]);

// 4. Build the Professional Bar Chart

var compChart = [Link](

accuracyComparison,

'Map Version',

['Overall Accuracy (%)']

.setChartType('ColumnChart')

.setOptions({

title: 'Overall Accuracy Comparison: Google Dynamic World vs. Custom RF Model',

hAxis: {title: 'LULC Map Version', textStyle: {bold: true}},

vAxis: {

title: 'Overall Accuracy (%)',

minValue: 0,

maxValue: 100,

textStyle: {bold: true}

},

// Make Dynamic World bars blue, and your Custom Model bar a bright contrasting orange/gold!

colors: ['#2ca02c'],

legend: {position: 'none'}

});

print('--- ACCURACY COMPARISON GRAPH ---');

print(compChart);

// ==============================================================================
// GRAPH: KAPPA COEFFICIENT COMPARISON (Dynamic World vs. My Custom Model)

// ==============================================================================

// 1. Get the Kappa coefficient of your Custom Model (Calculated in Step 2B)

var myKappa = [Link]();

// 2. Create a quick helper function to extract the Kappa for Dynamic World

var getDwKappa = function(dwImage) {

var dwTested = [Link]({

collection: allPoints, // Your hand-drawn ground truth pins

properties: ['class'],

scale: 10

});

return [Link]('class', 'dw_reference_class').kappa();

};

// Extract the exact Kappa values for the historical maps

var kap16 = getDwKappa(dw2016);

var kap20 = getDwKappa(dw2020);

var kap25 = getDwKappa(dw2025);

// 3. Bundle them together into a dataset for the Chart

var kappaComparison = [Link]([

[Link](null, {'Map Version': 'DW 2016', 'Kappa Coefficient': kap16}),

[Link](null, {'Map Version': 'DW 2020', 'Kappa Coefficient': kap20}),

[Link](null, {'Map Version': 'DW 2025', 'Kappa Coefficient': kap25}),

[Link](null, {'Map Version': 'My Custom Map', 'Kappa Coefficient': myKappa})

]);

// 4. Build the Professional Bar Chart

var kappaChart = [Link](

 kappaComparison,

'Map Version',

['Kappa Coefficient']

.setChartType('ColumnChart')

.setOptions({

title: 'Kappa Coefficient Comparison: Google DW vs. Custom RF Model',

hAxis: {title: 'LULC Map Version', textStyle: {bold: true}},

vAxis: {

title: 'Kappa Coefficient (0 to 1.0)',

minValue: 0,

maxValue: 1.0, // Kappa maxes out at 1.0!

textStyle: {bold: true}

},

// Using a rich purple color to differentiate it from the green Accuracy graph

colors: ['#9467bd'],

legend: {position: 'none'}

});

print('--- KAPPA COMPARISON GRAPH ---');

print(kappaChart);

// ==============================================================================

// STEP 4: Extract Land Surface Temperature (LST) for Multiple Years

// ==============================================================================

var getLST = function(year) {

var startDate = year + '-01-01';

var endDate = year + '-05-31';

var l8 = [Link]('LANDSAT/LC08/C02/T1_L2')

.filterBounds(roi)

.filterDate(startDate, endDate)
 .filter([Link]('CLOUD_COVER', 15))

.median()

.clip(roi);

var surfaceTemp = [Link]('ST_B10');

return [Link](0.00341802).add(149.0).subtract(273.15).rename('LST_Celsius_' +
year);

};

var lst2016 = getLST('2016');

var lst2020 = getLST('2020');

var lst2025 = getLST('2025'); // 2025 is now your master LST dataset

// ==============================================================================

// STEP 5: Calculate Average Temperature by LULC Class (Using 2025)

// ==============================================================================

var lstAndLulc = [Link](myCustomMap); // Using 2025 heat baseline!

var meanTempStats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'LULC_Class'}),

geometry: roi,

scale: 30,

maxPixels: 1e9

});

print('--- TEMPERATURE DATA (2025) ---');

print('Average LST by Land Cover Class (2025):', meanTempStats);

// ==============================================================================

// STEP 6: Visualize Everything on the Map

// ==============================================================================

[Link](landsat, {bands: ['SR_B4', 'SR_B3', 'SR_B2'], min: 0.0, max: 0.3}, 'Landsat True Color
(2025)');
var lulcPalette = ['red', 'green', 'blue', 'pink', 'yellow'];

[Link](myCustomMap, {min: 1, max: 5, palette: lulcPalette}, 'My LULC Map (2025)');

[Link](dw2016, {min: 1, max: 5, palette: lulcPalette}, 'Dynamic World (2016)', false);

[Link](dw2020, {min: 1, max: 5, palette: lulcPalette}, 'Dynamic World (2020)', false);

[Link](dw2025, {min: 1, max: 5, palette: lulcPalette}, 'Dynamic World (2025)', false);

var thermalPalette = ['blue', 'cyan', 'green', 'yellow', 'orange', 'red'];

[Link](lst2016, {min: 25, max: 40, palette: thermalPalette}, 'LST Map (2016)', false);

[Link](lst2020, {min: 25, max: 40, palette: thermalPalette}, 'LST Map (2020)', false);

[Link](lst2025, {min: 25, max: 40, palette: thermalPalette}, 'LST Map (2025)', false);

// ==============================================================================

// STEP 8A: Isolated Class Layers (Added to Default Layers Menu)

// ==============================================================================

// 1. Define our classes, values, and colors

var classList = [

{name: 'Urban', val: 1, color: 'red'},

{name: 'Forest', val: 2, color: 'green'},

{name: 'Water', val: 3, color: 'blue'},

{name: 'Vegetation', val: 4, color: 'pink'},

{name: 'Agriculture', val: 5, color: 'yellow'}

];

// 2. Define the three Dynamic World maps we want to slice up

var dwMaps = {'2016': dw2016, '2020': dw2020, '2025': dw2025};

// 3. Loop through every year and every class to create isolated layers

[Link](dwMaps).forEach(function(year) {

var dwMap = dwMaps[year];

 [Link](function(cls) {

// Isolate just this ONE class (e.g., only grab the 'Urban' pixels)

var isolatedClassMap = [Link]([Link]([Link]));

var layerName = year + ' - ' + [Link];

// Add it to the standard Layers dropdown, but keep it turned OFF (false) by default

[Link](isolatedClassMap, {min: [Link], max: [Link], palette: [[Link]]}, layerName, false);

});

});

// ==============================================================================

// STEP 8B: Professional UI Legends (LULC and LST)

// ==============================================================================

// ---------------------------------------------------------

// LEGEND 1: Land Cover Classes (Bottom-Left)

// ---------------------------------------------------------

var lulcLegend = [Link]({style: {position: 'bottom-left', padding: '10px 15px', backgroundColor:

'rgba(255, 255, 255, 0.9)'}});

var lulcTitle = [Link]({value: 'Land Cover Classes', style: {fontWeight: 'bold', fontSize: '15px', margin:
'0 0 6px 0'}});

[Link](lulcTitle);

var lulcColors = ['red', 'green', 'blue', 'pink', 'yellow'];

var lulcNames = ['Urban', 'Forest', 'Water', 'Vegetation', 'Agriculture'];

var makeLulcRow = function(color, name) {

var colorBox = [Link]({style: {backgroundColor: color, padding: '8px', margin: '0 8px 4px 0', border:
'1px solid black'}});

var description = [Link]({value: name, style: {margin: '0 0 4px 0'}});

return [Link]({widgets: [colorBox, description], layout: [Link]('horizontal')});

};

for (var i = 0; i < 5; i++) {

[Link](makeLulcRow(lulcColors[i], lulcNames[i]));

[Link](lulcLegend);

// ---------------------------------------------------------

// LEGEND 2: Surface Temperature Gradient (Bottom-Right)

// ---------------------------------------------------------

var lstLegend = [Link]({style: {position: 'bottom-right', padding: '10px 15px', backgroundColor:

'rgba(255, 255, 255, 0.9)', width: '250px'}});

var lstTitle = [Link]({value: 'Surface Temperature (°C)', style: {fontWeight: 'bold', fontSize: '15px',
margin: '0 0 6px 0'}});

[Link](lstTitle);

// Create the smooth color gradient using an Earth Engine thumbnail

var thermalPalette = ['blue', 'cyan', 'green', 'yellow', 'orange', 'red'];

var colorBar = [Link]({

image: [Link]().select(0),

params: {

bbox: [0, 0, 1, 0.1], // Creates a horizontal bar

dimensions: '200x15', // Width and height of the bar

format: 'png',

min: 0,

max: 1,

palette: thermalPalette,

},

style: {stretch: 'horizontal', margin: '0px 0px', border: '1px solid black'}

});

[Link](colorBar);
// Add the temperature numbers below the gradient bar

var legendLabels = [Link]({

widgets: [

[Link]('25°C', {margin: '4px 0px 0px 0px', fontSize: '13px', fontWeight: 'bold'}),

[Link]('32°C', {margin: '4px 0px 0px 0px', fontSize: '13px', fontWeight: 'bold', textAlign: 'center',
stretch: 'horizontal'}),

[Link]('40°C+', {margin: '4px 0px 0px 0px', fontSize: '13px', fontWeight: 'bold'})

],

layout: [Link]('horizontal')

});

[Link](legendLabels);

[Link](lstLegend);

// ==============================================================================

// STEP 9: Cloud Cover Data Panel

// ==============================================================================

var cloudPanel = [Link]({

style: {position: 'top-right', padding: '8px 15px', backgroundColor: 'white'}

});

var cloudLabel = [Link]({

value: 'Calculating Cloud Cover...',

style: {fontWeight: 'bold', fontSize: '14px', margin: '0'}

});

[Link](cloudLabel);

[Link](cloudPanel);

var avgCloudCover = landsatCollection.aggregate_mean('CLOUD_COVER');

[Link](function(cloudPercentage){

if (cloudPercentage !== null) {

 [Link]('Average Cloud Cover (2025): ' + [Link](2) + '%');

} else {

[Link]('No Cloud Data Found');

});

// ==============================================================================

// STEP 9: VALIDATE THE TEMPERATURE MODEL (RMSE & R-Squared)

// ==============================================================================

print('=========================================');

print('--- TEMPERATURE MODEL ACCURACY ---');

// 0. REBUILD THE LST TREND (To prevent 'undefined' errors)

var buildTrendData = function(year) {

var startDate = [Link](year, 1, 1);

var endDate = [Link](year, 5, 31);

var l8 = [Link]('LANDSAT/LC08/C02/T1_L2')

.filterBounds(roi)

.filterDate(startDate, endDate)

.filter([Link]('CLOUD_COVER', 15))

.median();

var lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15);

return [Link](year).rename('year').addBands([Link]('temp')).toFloat();

};

var years = [Link](2014, 2025);

var trendCollection = [Link]([Link](buildTrendData));

var linearFit = [Link]([Link]());

// 1. Get ACTUAL 2025 Temperature

var actualLST2025 = buildTrendData(2025).select('temp').rename('LST').clip(roi);

// 2. Predict 2025 Temperature using the model

var predictedLST2025 = [Link]('scale').multiply(2025)

.add([Link]('offset'))

.rename('LST')

.clip(roi);

// --- A. CALCULATE RMSE ---

var error = [Link](actualLST2025);

var squaredError = [Link](2);

var mseDictionary = [Link]({

reducer: [Link](),

geometry: roi,

scale: 30,

maxPixels: 1e9

});

var rmse = [Link]([Link]('LST')).sqrt();

print('RMSE (°C):', rmse);

// --- B. CALCULATE R-SQUARED (The "Accuracy %") ---

// Combine the predicted and actual maps into one image for comparison

var combinedLST =
[Link]('predicted').addBands([Link]('actual'));

// Calculate the Pearson's Correlation Coefficient (r)

var correlationDict = [Link]({

reducer: [Link](),

geometry: roi,

scale: 30,
 maxPixels: 1e9

});

// Extract 'correlation', square it, and convert to percentage!

var r = [Link]([Link]('correlation'));

var rSquared = [Link](2);

var accuracyPercent = [Link](100);

print('Pearson Correlation (r):', r);

print('R-Squared (R²):', rSquared);

print('Temperature Accuracy %:', accuracyPercent);

print('=========================================');

// ==============================================================================

// STEP 10: Multi-Year Temperature Chart for Active Classes (2016, 2020, 2025)

// ==============================================================================

var tempYears = [Link]([2016, 2020, 2025]);

var getTempByYearChart = function(year) {

var y = [Link](year);

var startDate = [Link](y, 1, 1);

var endDate = [Link](y, 5, 31);

var l8 = [Link]('LANDSAT/LC08/C02/T1_L2')

.filterBounds(roi)

.filterDate(startDate, endDate)

.filter([Link]('CLOUD_COVER', 15))

.median();

var lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15);

 var urbanTemp = [Link]([Link](1)).reduceRegion([Link](), roi, 30,
null, null, false, 1e9).get('ST_B10');

var forestTemp = [Link]([Link](2)).reduceRegion([Link](), roi, 30,

null, null, false, 1e9).get('ST_B10');

var waterTemp = [Link]([Link](3)).reduceRegion([Link](), roi, 30,

null, null, false, 1e9).get('ST_B10');

var vegTemp = [Link]([Link](4)).reduceRegion([Link](), roi, 30, null,

null, false, 1e9).get('ST_B10');

var agriTemp = [Link]([Link](5)).reduceRegion([Link](), roi, 30,

null, null, false, 1e9).get('ST_B10');

return [Link](null, {

'Year': [Link]([Link]()),

'Urban': urbanTemp,

'Forest': forestTemp,

'Water': waterTemp,

'Vegetation': vegTemp,

'Agriculture': agriTemp

});

};

var yearlyTemps = [Link]([Link](getTempByYearChart));

var tempChart = [Link](

yearlyTemps,

'Year',

['Urban', 'Forest', 'Water', 'Vegetation', 'Agriculture']

).setOptions({

title: 'Land Surface Temperature Trend by Class (2016, 2020, 2025)',

hAxis: {title: 'Year'},

vAxis: {title: 'Average Surface Temperature (°C)'},

colors: ['red', 'green', 'blue', 'pink', 'yellow'],

lineWidth: 3,
 pointSize: 7

});

print('--- MULTI-YEAR TEMPERATURE TREND ---');

print(tempChart);

// ==============================================================================

// STEP 12: Multi-Year Land Cover Area Chart (2016, 2020, 2025)

// ==============================================================================

var dwYears = [Link]([2016, 2020, 2025]);

var getAreaByYear = function(year) {

var y = [Link](year);

var dwMap = getDynamicWorld(y);

var pixelAreaSqKm = [Link]().divide(1e6);

var urbanArea = [Link]([Link](1)).reduceRegion([Link](), roi,

10, null, null, false, 1e9).get('area');

var forestArea = [Link]([Link](2)).reduceRegion([Link](), roi,

10, null, null, false, 1e9).get('area');

var waterArea = [Link]([Link](3)).reduceRegion([Link](), roi,

10, null, null, false, 1e9).get('area');

var vegArea = [Link]([Link](4)).reduceRegion([Link](), roi, 10,

null, null, false, 1e9).get('area');

var agriArea = [Link]([Link](5)).reduceRegion([Link](), roi, 10,

null, null, false, 1e9).get('area');

return [Link](null, {

'Year': [Link]([Link]()),

'Urban Area (Sq Km)': urbanArea,

'Forest Area (Sq Km)': forestArea,

'Water Area (Sq Km)': waterArea,

'Vegetation Area (Sq Km)': vegArea,

 'Agriculture Area (Sq Km)': agriArea

});

};

var yearlyAreas = [Link]([Link](getAreaByYear));

var areaChart = [Link](

yearlyAreas,

'Year',

['Urban Area (Sq Km)', 'Forest Area (Sq Km)', 'Water Area (Sq Km)', 'Vegetation Area (Sq Km)',
'Agriculture Area (Sq Km)']

).setChartType('ColumnChart')

.setOptions({

title: 'Urban Expansion vs Natural Land Loss (2016, 2020, 2025)',

hAxis: {title: 'Year'},

vAxis: {title: 'Total Area (Square Kilometers)'},

colors: ['red', 'green', 'blue', 'pink', 'yellow']

});

print('--- MULTI-YEAR LAND COVER AREA ---');

print(areaChart);

// ==============================================================================

// STEP 13: Predict Land Cover Area for 2050 (2016 to 2025 Trendline)

// ==============================================================================

// STEP 14: Predict 2050 Land Surface Temperature (Linear Regression)

// ==============================================================================

// ==============================================================================
// STEP 15: Multi-Decade Spatial Prediction for 2030, 2040, and 2050

// ==============================================================================

print('Calculating Multi-Decade Urban Sprawl Trend... (This may take a minute)');

var classList = [Link]([1, 2, 3, 4, 5]);

var yearsListDW = [Link](2016, 2025);

// The kernel radius determines how far the "sprawl" virus can look around itself

var neighborhood = [Link]({radius: 500, units: 'meters'});

// 1. Get the Linear Fit (Slope and Offset) for ONE class over 10 years

var getLinearFitForClass = function(classNum) {

var cNum = [Link](classNum);

var yearlyProbabilities = [Link]([Link](function(year) {

var y = [Link](year);

var dwMap = getDynamicWorld(y);

var binaryMap = [Link](cNum);

var spatialDensity = [Link]({

reducer: [Link](),

kernel: neighborhood

});

return [Link](y).rename('year')

.addBands([Link]('density')).toFloat();

}));

// Return the mathematical slope (scale) and intercept (offset) for this class

return [Link]([Link]());

};
// 2. Calculate the regression coefficients for ALL 5 classes ONCE to save memory

var linearFits = [Link](getLinearFitForClass);

// 3. Master function to predict ANY future year based on the saved slope

var predictFutureLULC = function(targetYear) {

var y = [Link](targetYear);

var futureProbabilities = [Link]([Link](function(classNum) {

var cNum = [Link](classNum);

var fit = [Link]([Link]([Link](1)));

// y = mx + b (Slope * Year + Intercept)

var predictedDensity = [Link]('scale').multiply(y)

.add([Link]('offset'))

.max(0);

return [Link]([Link]('c').cat([Link]('%d')));

})).toBands();

// Find the winning class for every pixel

return [Link]().arrayArgmax().arrayGet([0]).add(1)

.rename('Predicted_LULC_' + targetYear)

.clip(roi);

};

// 4. Generate the 2030, 2040, and 2050 Maps!

var predictedMap2030 = predictFutureLULC(2030);

var predictedMap2040 = predictFutureLULC(2040);

var predictedMap2050_MultiYear = predictFutureLULC(2050);

// 5. Add them to the map viewer

[Link](predictedMap2030, {min: 1, max: 5, palette: lulcPalette}, 'Predicted Map (2030)',

false);

[Link](predictedMap2040, {min: 1, max: 5, palette: lulcPalette}, 'Predicted Map (2040)',

false);

[Link](predictedMap2050_MultiYear, {min: 1, max: 5, palette: lulcPalette}, 'Predicted Map

(2050)', false);

print('Multi-Decade Spatial Prediction Complete!');

// ==============================================================================

// EXPORT TRAINING POINTS TO ASSETS (For Google Colab)

// ==============================================================================

// 1. Combine all your manual pins into one master file

var allPoints = [Link](forest)

.merge(water)

.merge(vegetation)

.merge(agriculture);

// 2. Export that master file directly into your Google Earth Engine Assets

[Link]({

collection: allPoints,

description: 'Export_My_Training_Points',

assetId: 'projects/fourth-vehicle-452506-f2/assets/training_points' // Saving right next to your

shape1!

});

// ==============================================================================

// FEATURE EXTRACTED FROM PAPER: 3-Part Seasonal NDBI vs LST (Figure 9)

// ==============================================================================

print('Generating Seasonal NDBI vs LST Charts (Jan, Jun, Oct)...');

// 1. Create a function that isolates a SPECIFIC SEASON across the three years
// 1. Create a function that isolates a SPECIFIC SEASON across the three years

var getSeasonalBinnedNDBIData = function(year, monthNumber) {

var startDate = [Link](year, monthNumber, 1);

// FIX 1: Expand to a full 3-month season to survive the monsoon!

var endDate = [Link](3, 'month');

// FIX 2: Merge the High Quality (Tier 1) and Cloudy Backup (Tier 2) collections

var t1 = [Link]('LANDSAT/LC08/C02/T1_L2');

var t2 = [Link]('LANDSAT/LC08/C02/T2_L2');

var l8 = [Link](t2)

.filterBounds(roi)

.filterDate(startDate, endDate)

.median()

.clip(roi);

// Calculate LST

var lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

// Calculate NDBI and Bin it

var ndbi = [Link](['SR_B6', 'SR_B5']);

var ndbiBin = [Link](10).round().divide(10).rename('NDBI_Bin');

var combined = [Link](ndbiBin);

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'NDBI_Bin'}),

geometry: roi,

scale: 100,

maxPixels: 1e10
 });

var groups = [Link]([Link]('groups'));

var features = [Link](function(g) {

var dict = [Link](g);

return [Link](null, {

'NDBI_Bin': [Link]('NDBI_Bin'),

'Mean_LST': [Link]('mean'),

'Year': [Link]([Link]())

});

});

return [Link](features);

};

// 2. Build a Helper Function to generate and print the NDBI Chart for a given month

var printSeasonalNDBIChart = function(monthNum, monthName) {

var data16 = getSeasonalBinnedNDBIData(2016, monthNum);

var data20 = getSeasonalBinnedNDBIData(2020, monthNum);

var data25 = getSeasonalBinnedNDBIData(2025, monthNum);

var allData = [Link](data20).merge(data25);

var chart = [Link]({

features: allData,

xProperty: 'NDBI_Bin',

yProperty: 'Mean_LST',

seriesProperty: 'Year'

})

.setChartType('LineChart')

.setOptions({

title: 'Relationship between Mean LST and NDBI (' + monthName + ')',
 hAxis: {title: 'NDBI Value (Higher = More Concrete/Built-up)', minValue: -0.6, maxValue: 1.0},

vAxis: {title: 'Mean LST (°C)'},

lineWidth: 3,

pointSize: 4,

series: {

0: {color: '#2ca02c', lineDashStyle: [4, 4]}, // 2016

1: {color: '#1f77b4', lineDashStyle: [2, 2]}, // 2020

2: {color: '#d62728', lineDashStyle: [10, 2]} // 2025

},

curveType: 'function'

});

print('--- NDBI vs LST (' + [Link]() + ') ---');

print(chart);

};

// 3. Command Earth Engine to print all three seasonal charts!

printSeasonalNDBIChart(1, 'January');

printSeasonalNDBIChart(6, 'June');

printSeasonalNDBIChart(10, 'October');

// ==============================================================================

// FEATURE EXTRACTED FROM PAPER: 3-Part Seasonal NDVI vs LST (Figure 7)

// ==============================================================================

print('Generating Seasonal NDVI vs LST Charts (Jan, Jun, Oct)...');

// 1. Create a function that isolates a SPECIFIC SEASON across the three years for Vegetation (NDVI)

var getSeasonalBinnedNDVIData = function(year, monthNumber) {

var startDate = [Link](year, monthNumber, 1);

// Using the 3-month window to survive the June monsoon!

var endDate = [Link](3, 'month');

// Merge the High Quality (Tier 1) and Cloudy Backup (Tier 2) collections

var t1 = [Link]('LANDSAT/LC08/C02/T1_L2');

var t2 = [Link]('LANDSAT/LC08/C02/T2_L2');

var l8 = [Link](t2)

.filterBounds(roi)

.filterDate(startDate, endDate)

.median()

.clip(roi);

// Calculate LST

var lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

// Calculate NDVI (Normalized Difference Vegetation Index)

// Formula for Landsat 8: (NIR - Red) / (NIR + Red) -> (Band 5 and Band 4)

var ndvi = [Link](['SR_B5', 'SR_B4']);

// Bin the NDVI values (rounding them into clean categories like 0.1, 0.2, 0.3)

var ndviBin = [Link](10).round().divide(10).rename('NDVI_Bin');

// Earth Engine strictly requires the grouping band to be added LAST (Band 1).

var combined = [Link](ndviBin);

// Calculate the average temperature for each vegetation bin

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'NDVI_Bin'}),

geometry: roi,

scale: 100,
 maxPixels: 1e10

});

// Format for the chart

var groups = [Link]([Link]('groups'));

var features = [Link](function(g) {

var dict = [Link](g);

return [Link](null, {

'NDVI_Bin': [Link]('NDVI_Bin'),

'Mean_LST': [Link]('mean'),

'Year': [Link]([Link]())

});

});

return [Link](features);

};

// 2. Build a Helper Function to generate and print the NDVI Chart for a given month

var printSeasonalNDVIChart = function(monthNum, monthName) {

var data16 = getSeasonalBinnedNDVIData(2016, monthNum);

var data20 = getSeasonalBinnedNDVIData(2020, monthNum);

var data25 = getSeasonalBinnedNDVIData(2025, monthNum);

var allData = [Link](data20).merge(data25);

var chart = [Link]({

features: allData,

xProperty: 'NDVI_Bin',

yProperty: 'Mean_LST',

seriesProperty: 'Year'

})

.setChartType('LineChart')
 .setOptions({

title: 'Relationship between Mean LST and NDVI (' + monthName + ')',

hAxis: {title: 'NDVI Value (Higher = Denser Vegetation)', minValue: -0.6, maxValue: 1.0},

vAxis: {title: 'Mean LST (°C)'},

lineWidth: 3,

pointSize: 4,

// Mimicking the paper's exact line styles

series: {

0: {color: '#2ca02c', lineDashStyle: [4, 4]}, // 2016 (Green Dash)

1: {color: '#1f77b4', lineDashStyle: [2, 2]}, // 2020 (Blue Dot)

2: {color: '#d62728', lineDashStyle: [10, 2]} // 2025 (Red Long Dash)

},

curveType: 'function'

});

print('--- NDVI vs LST (' + [Link]() + ') ---');

print(chart);

};

// 3. Command Earth Engine to print all three seasonal charts!

printSeasonalNDVIChart(1, 'January');

printSeasonalNDVIChart(6, 'June');

printSeasonalNDVIChart(10, 'October');

// ==============================================================================

// FEATURE EXTRACTED FROM PAPER: 3-Part Seasonal MNDWI vs LST (Figure 10)

// ==============================================================================

print('Generating Seasonal MNDWI vs LST Charts (Jan, Jun, Oct)...');

// 1. Create a function that isolates a SPECIFIC SEASON across the three years for Water (MNDWI)

var getSeasonalBinnedMNDWIData = function(year, monthNumber) {

var startDate = [Link](year, monthNumber, 1);

// Using the 3-month window to survive the monsoon!

var endDate = [Link](3, 'month');

// Merge the High Quality (Tier 1) and Cloudy Backup (Tier 2) collections

var t1 = [Link]('LANDSAT/LC08/C02/T1_L2');

var t2 = [Link]('LANDSAT/LC08/C02/T2_L2');

var l8 = [Link](t2)

.filterBounds(roi)

.filterDate(startDate, endDate)

.median()

.clip(roi);

// Calculate LST

var lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

// Calculate MNDWI (Modified Normalized Difference Water Index)

// NOTE: To match the paper's visuals/text where -1 is water and +1 is dry land,

// we use (SWIR1 - Green) instead of the standard (Green - SWIR1).

// Landsat 8: Band 6 (SWIR1) and Band 3 (Green)

var mndwi = [Link](['SR_B6', 'SR_B3']);

// Bin the MNDWI values (rounding them into clean categories like 0.1, 0.2, 0.3)

var mndwiBin = [Link](10).round().divide(10).rename('MNDWI_Bin');

// Earth Engine strictly requires the grouping band to be added LAST (Band 1).

var combined = [Link](mndwiBin);

// Calculate the average temperature for each water bin

 var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'MNDWI_Bin'}),

geometry: roi,

scale: 100,

maxPixels: 1e10

});

// Format for the chart

var groups = [Link]([Link]('groups'));

var features = [Link](function(g) {

var dict = [Link](g);

return [Link](null, {

'MNDWI_Bin': [Link]('MNDWI_Bin'),

'Mean_LST': [Link]('mean'),

'Year': [Link]([Link]())

});

});

return [Link](features);

};

// 2. Build a Helper Function to generate and print the MNDWI Chart for a given month

var printSeasonalMNDWIChart = function(monthNum, monthName) {

var data16 = getSeasonalBinnedMNDWIData(2016, monthNum);

var data20 = getSeasonalBinnedMNDWIData(2020, monthNum);

var data25 = getSeasonalBinnedMNDWIData(2025, monthNum);

var allData = [Link](data20).merge(data25);

var chart = [Link]({

features: allData,

xProperty: 'MNDWI_Bin',
 yProperty: 'Mean_LST',

seriesProperty: 'Year'

})

.setChartType('LineChart')

.setOptions({

title: 'Relationship between Mean LST and MNDWI (' + monthName + ')',

hAxis: {title: 'MNDWI Value (Higher = Drier/Impervious, Lower = Water)', minValue: -0.8,
maxValue: 1.0},

vAxis: {title: 'Mean LST (°C)'},

lineWidth: 3,

pointSize: 4,

// Mimicking the paper's exact line styles

series: {

0: {color: '#2ca02c', lineDashStyle: [4, 4]}, // 2016 (Green Dash)

1: {color: '#1f77b4', lineDashStyle: [2, 2]}, // 2020 (Blue Dot)

2: {color: '#d62728', lineDashStyle: [10, 2]} // 2025 (Red Long Dash)

},

curveType: 'function'

});

print('--- MNDWI vs LST (' + [Link]() + ') ---');

print(chart);

};

// 3. Command Earth Engine to print all three seasonal charts!

printSeasonalMNDWIChart(1, 'January');

printSeasonalMNDWIChart(6, 'June');

printSeasonalMNDWIChart(10, 'October');

// ==============================================================================

// FEATURE EXTRACTED FROM PAPER: 3-Part Seasonal NDMI vs LST (Figure 8)

// ==============================================================================

print('Generating Seasonal NDMI vs LST Charts (Jan, Jun, Oct)...');

// 1. Create a function that isolates a SPECIFIC SEASON across the three years for Moisture (NDMI)

var getSeasonalBinnedNDMIData = function(year, monthNumber) {

var startDate = [Link](year, monthNumber, 1);

// Shift the summer window to catch pre-monsoon heat

if (monthNumber === 6) {

startDate = [Link](year, 5, 15);

var endDate = [Link](3, 'month');

// Dual-Satellite Fusion (Landsat 8 + Landsat 9) for maximum cloud-free chances

var l8_t1 = [Link]('LANDSAT/LC08/C02/T1_L2');

var l8_t2 = [Link]('LANDSAT/LC08/C02/T2_L2');

var l9_t1 = [Link]('LANDSAT/LC09/C02/T1_L2');

var l9_t2 = [Link]('LANDSAT/LC09/C02/T2_L2');

var combinedLandsat = l8_t1.merge(l8_t2).merge(l9_t1).merge(l9_t2);

var lSat = combinedLandsat

.filterBounds(roi)

.filterDate(startDate, endDate)

.map(function(img) {

var qa = [Link]('QA_PIXEL');

var cloudBitMask = 1 << 3;

var shadowBitMask = 1 << 4;

var cleanPixels = [Link](cloudBitMask).eq(0)

.and([Link](shadowBitMask).eq(0));
 var realDataMask = [Link]('ST_B10').gt(0);

return [Link](cleanPixels).updateMask(realDataMask);

})

.median()

.clip(roi);

// Calculate LST

var lst =
[Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

// Calculate NDMI (Normalized Difference Moisture Index)

// Landsat 8/9 Formula: (NIR - SWIR1) / (NIR + SWIR1) -> (Band 5 and Band 6)

var ndmi = [Link](['SR_B5', 'SR_B6']);

var ndmiBin = [Link](10).round().divide(10).rename('NDMI_Bin');

var combined = [Link](ndmiBin);

// ABSOLUTE SAFETY NET: Throw out any broken pixels colder than 15°C!

combined = [Link]([Link]('Mean_LST').gt(15));

// Calculate the average temperature for each moisture bin

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'NDMI_Bin'}),

geometry: roi,

scale: 100,

maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

var features = [Link](function(g) {

var dict = [Link](g);

 return [Link](null, {

'NDMI_Bin': [Link]('NDMI_Bin'),

'Mean_LST': [Link]('mean'),

'Year': [Link]([Link]())

});

});

return [Link](features);

};

// 2. Build a Helper Function to generate and print the NDMI Chart

var printSeasonalNDMIChart = function(monthNum, monthName) {

var data16 = getSeasonalBinnedNDMIData(2016, monthNum);

var data20 = getSeasonalBinnedNDMIData(2020, monthNum);

var data25 = getSeasonalBinnedNDMIData(2025, monthNum);

var allData = [Link](data20).merge(data25);

var chart = [Link]({

features: allData,

xProperty: 'NDMI_Bin',

yProperty: 'Mean_LST',

seriesProperty: 'Year'

})

.setChartType('LineChart')

.setOptions({

title: 'Relationship between Mean LST and NDMI (' + monthName + ')',

hAxis: {title: 'NDMI Value (Higher = More Moisture)', minValue: -0.8, maxValue: 0.8},

vAxis: {title: 'Mean LST (°C)'},

lineWidth: 3,

pointSize: 4,

series: {
 0: {color: '#2ca02c', lineDashStyle: [4, 4]}, // 2016

1: {color: '#1f77b4', lineDashStyle: [2, 2]}, // 2020

2: {color: '#d62728', lineDashStyle: [10, 2]} // 2025

},

curveType: 'function'

});

print('--- NDMI vs LST (' + [Link]() + ') ---');

print(chart);

};

// 3. Command Earth Engine to print all three seasonal charts!

printSeasonalNDMIChart(1, 'January');

printSeasonalNDMIChart(6, 'June');

printSeasonalNDMIChart(10, 'October');

// ==============================================================================

// FEATURE EXTRACTED FROM PAPER: Mean LST by LULC Class (Figure 5)

// ==============================================================================

print('Generating Seasonal LST by Class Charts (Corrected for Clouds)...');

// 1. Helper function to safely extract LST with a Cloud/Zero Mask

// 1. Helper function to safely extract LST with a Cloud/Zero Mask

// 1. Helper function to safely extract LST with a Strict QA Cloud Mask

// 1. Helper function: Dual-Satellite (Landsat 8 + Landsat 9) with QA Cloud Mask

var getMonthlyLST = function(year, monthNumber) {

var startDate = [Link](year, monthNumber, 1);

// FIX: Shift June window to pre-monsoon (Mid-May) to capture peak heat

// June itself is heavily cloud-covered due to monsoon onset

if (monthNumber === 6) {
 startDate = [Link](year, 5, 1); // Start from May 1st

var endDate = [Link](3, 'month');

// ... rest of function unchanged

// Landsat 8 (Tier 1 & Tier 2)

var l8_t1 = [Link]('LANDSAT/LC08/C02/T1_L2');

var l8_t2 = [Link]('LANDSAT/LC08/C02/T2_L2');

// FIX: Add Landsat 9 (Tier 1 & Tier 2) to DOUBLE our observation frequency for 2025!

var l9_t1 = [Link]('LANDSAT/LC09/C02/T1_L2');

var l9_t2 = [Link]('LANDSAT/LC09/C02/T2_L2');

// Merge all 4 collections into one massive dataset

var combinedLandsat = l8_t1.merge(l8_t2).merge(l9_t1).merge(l9_t2);

var lSat = combinedLandsat

.filterBounds(roi)

.filterDate(startDate, endDate)

// Apply the strict NASA QA_PIXEL mask

.map(function(img) {

var qa = [Link]('QA_PIXEL');

// Bit 3 is Clouds, Bit 4 is Cloud Shadows

var cloudBitMask = 1 << 3;

var shadowBitMask = 1 << 4;

var cleanPixels = [Link](cloudBitMask).eq(0)

.and([Link](shadowBitMask).eq(0));

var realDataMask = [Link]('ST_B10').gt(0);

return [Link](cleanPixels).updateMask(realDataMask);

})
 .median()

.clip(roi);

// Convert to Celsius

var lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15);

// Final safety net: Mumbai surface temp rarely drops below 22°C during these months

var finalMask = [Link](22);

return [Link](finalMask);

};

// 2. Master function to generate the 3-Month Bar Chart

var generateFigure5 = function(chartTitle, lstYear, lulcMap) {

var lstJan = getMonthlyLST(lstYear, 1).rename('January');

var lstJun = getMonthlyLST(lstYear, 6).rename('June');

var lstOct = getMonthlyLST(lstYear, 10).rename('October');

var seasonalLST = [Link](lstJun).addBands(lstOct);

var lulcForChart = [Link](1).rename('class');

// 3. Build the Professional Grouped Column Chart

var chart = [Link]({

image: [Link](lulcForChart),

classBand: 'class',

region: roi,

reducer: [Link](),

scale: 150,

classLabels: ['Urban', 'Forest', 'Water', 'Vegetation', 'Agriculture']

})

.setChartType('ColumnChart')
 .setOptions({

title: 'Mean LST by LULC Class: ' + chartTitle,

hAxis: {title: 'Land Cover Class', textStyle: {bold: true}},

vAxis: {title: 'Mean LST (°C)', textStyle: {bold: true}},

colors: ['#1f77b4', '#d62728', '#ff7f0e']

});

print('--- ' + [Link]() + ' ---');

print(chart);

};

// 3. Command Earth Engine to generate the 4 charts

generateFigure5('Dynamic World 2016', 2016, dw2016);

generateFigure5('Dynamic World 2020', 2020, dw2020);

generateFigure5('Dynamic World 2025', 2025, dw2025);

// ==============================================================================

// FEATURE EXTRACTED FROM PAPER: LULC Area & Change Statistics (Table 2)

// ==============================================================================

print('Calculating Land Transformation Statistics (This may take a moment)...');

// 1. Define your 5 classes

var classNames = [Link](['1. Urban', '2. Forest', '3. Water', '4. Vegetation', '5. Agriculture']);

var classValues = [Link]([1, 2, 3, 4, 5]);

// 2. Helper Function: Calculate Area (in Hectares) for every class in a map

var getMapAreas = function(lulcImage) {

// Multiply pixel area (m2) by the map classes, then divide by 10,000 for Hectares

var areaImage = [Link]().divide(10000).addBands([Link]('class'));

var stats = [Link]({

 reducer: [Link]().group({groupField: 1, groupName: 'class'}),

geometry: roi,

scale: 30, // Using 30m resolution to match Landsat scale

maxPixels: 1e10

});

// Extract the results into a Dictionary: {'1': area, '2': area, etc.}

var groups = [Link]([Link]('groups'));

var classAreas = [Link](function(item) {

var d = [Link](item);

var classString = [Link]([Link]('class')).toInt().format('%d');

return [Link]([classString, [Link]('sum')]);

});

// Ensure all 5 classes exist in the dictionary, even if area is 0

var defaultDict = [Link]({'1':0, '2':0, '3':0, '4':0, '5':0});

return [Link]([Link]([Link]()), true);

};

// 3. Extract the Hectares for 2016, 2020, and 2025

var area16 = getMapAreas(dw2016);

var area20 = getMapAreas(dw2020);

var area25 = getMapAreas(dw2025); // Using your custom map for the present!

// 4. Calculate the Total City Area to figure out percentages

var totalArea = [Link]([Link]().reduce([Link]()));

// ==============================================================================

// FEATURE EXTRACTED FROM PAPER: LULC Area & Change Statistics (Table 1)

// ==============================================================================

print('Calculating Land Transformation Statistics (Sq Km)...');

// 1. Define your 5 classes

var classNames = [Link](['1. Urban', '2. Forest', '3. Water', '4. Vegetation', '5. Agriculture']);

var classValues = [Link]([1, 2, 3, 4, 5]);

// 2. Helper Function: Calculate Area (in Sq Km) for every class in a map

var getMapAreas = function(lulcImage) {

// FIX: Divide pixel area by 1,000,000 to convert square meters into Square Kilometers

var areaImage = [Link]().divide(1000000).addBands([Link]('class'));

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'class'}),

geometry: roi,

scale: 30,

maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

var classAreas = [Link](function(item) {

var d = [Link](item);

var classString = [Link]([Link]('class')).toInt().format('%d');

return [Link]([classString, [Link]('sum')]);

});

var defaultDict = [Link]({'1':0, '2':0, '3':0, '4':0, '5':0});

return [Link]([Link]([Link]()), true);

};

// 3. Extract the Sq Km for 2016, 2020, and 2025

// 3. Extract the Sq Km for 2016, 2020, and 2025

// 3. Extract the Sq Km for Past, Present, and FUTURE maps!

var area16 = getMapAreas(dw2016);

var area20 = getMapAreas(dw2020);

// The 2025 Correction (Once Urban, Always Urban)

var dw2025_corrected = [Link]([Link](1), 1);

var area25 = getMapAreas(dw2025_corrected);

// FUTURE PREDICTIONS: Apply the irreversibility filter so future concrete never vanishes!

// (Make sure these variable names match exactly what you named them in your prediction code)

var pred30 = [Link](dw2025_corrected.eq(1), 1);

var pred40 = [Link]([Link](1), 1);

var pred50 = predictedMap2050_MultiYear.where([Link](1), 1);

var area30 = getMapAreas(pred30);

var area40 = getMapAreas(pred40);

var area50 = getMapAreas(pred50);

// 4. Calculate the Total City Area

var totalArea = [Link]([Link]().reduce([Link]()));

// 5. Build the rows for our Mega-Table

var tableRows = [Link](function(c) {

var cNum = [Link](c);

var cKey = [Link]('%d');

var name = [Link]([Link](1));

// Get raw areas in sq km for ALL years

var a16 = [Link]([Link](cKey));

var a20 = [Link]([Link](cKey));

var a25 = [Link]([Link](cKey));

var a30 = [Link]([Link](cKey));

var a40 = [Link]([Link](cKey));

 var a50 = [Link]([Link](cKey));

// Calculate Future Change in Area (Sq Km)

var diff25_30 = [Link](a25);

var diff30_40 = [Link](a30);

var diff40_50 = [Link](a40);

// Calculate the MASSIVE total change from the beginning of your study to the end

var diff_Total_16_to_50 = [Link](a16);

// Create strict column names for the database

return [Link](null, {

'LULC_Class': name,

'Area_2016': [Link]('%.2f'),

'Area_2025': [Link]('%.2f'),

'Area_2030': [Link]('%.2f'),

'Area_2040': [Link]('%.2f'),

'Area_2050': [Link]('%.2f'),

'Change_25_to_30': diff25_30.format('%.2f').cat(' sq km'),

'Change_30_to_40': diff30_40.format('%.2f').cat(' sq km'),

'Change_40_to_50': diff40_50.format('%.2f').cat(' sq km'),

'TOTAL_CHANGE_16_to_50': diff_Total_16_to_50.format('%.2f').cat(' sq km')

});

});

// 6. Generate the Interactive Console Table

var statsTable = [Link]([Link](tableRows), 'LULC_Class')

.setChartType('Table')

.setOptions({

title: 'Multi-Decade LULC Transformation & Prediction (2016 - 2050)',

allowHtml: true,
 pageSize: 5

});

print('--- FUTURE PREDICTION AREA STATISTICS ---');

print(statsTable);

// ==============================================================================

// STEP 17: Multi-Decade Predicted LST Bar Chart (2030, 2040, 2050)

// ==============================================================================

print('Generating Future Temperature Bar Chart...');

// 1. Project the physical temperature maps for the specific future decades

// (This safely uses the 'linearFit' regression you calculated in Step 14)

var projected2030Map =
[Link]('scale').multiply(2030).add([Link]('offset')).rename('LST_2030').clip(roi);

var projected2040Map =
[Link]('scale').multiply(2040).add([Link]('offset')).rename('LST_2040').clip(roi);

var projected2050Map =
[Link]('scale').multiply(2050).add([Link]('offset')).rename('LST_2050').clip(roi);

// 2. Helper function to extract temperature using the matching FUTURE land cover map

var extractFutureTemp = function(tempImage, lulcMap, classNum, bandName) {

return [Link]([Link](classNum))

.reduceRegion({

reducer: [Link](),

geometry: roi,

scale: 30,

maxPixels: 1e9

}).get(bandName);

};
// 3. Extract the temperatures for all 5 classes across all 3 future decades

// URBAN

var u30_t = extractFutureTemp(projected2030Map, predictedMap2030, 1, 'LST_2030');

var u40_t = extractFutureTemp(projected2040Map, predictedMap2040, 1, 'LST_2040');

var u50_t = extractFutureTemp(projected2050Map, predictedMap2050_MultiYear, 1, 'LST_2050');

// FOREST

var f30_t = extractFutureTemp(projected2030Map, predictedMap2030, 2, 'LST_2030');

var f40_t = extractFutureTemp(projected2040Map, predictedMap2040, 2, 'LST_2040');

var f50_t = extractFutureTemp(projected2050Map, predictedMap2050_MultiYear, 2, 'LST_2050');

// WATER

var w30_t = extractFutureTemp(projected2030Map, predictedMap2030, 3, 'LST_2030');

var w40_t = extractFutureTemp(projected2040Map, predictedMap2040, 3, 'LST_2040');

var w50_t = extractFutureTemp(projected2050Map, predictedMap2050_MultiYear, 3, 'LST_2050');

// VEGETATION

var v30_t = extractFutureTemp(projected2030Map, predictedMap2030, 4, 'LST_2030');

var v40_t = extractFutureTemp(projected2040Map, predictedMap2040, 4, 'LST_2040');

var v50_t = extractFutureTemp(projected2050Map, predictedMap2050_MultiYear, 4, 'LST_2050');

// AGRICULTURE

var a30_t = extractFutureTemp(projected2030Map, predictedMap2030, 5, 'LST_2030');

var a40_t = extractFutureTemp(projected2040Map, predictedMap2040, 5, 'LST_2040');

var a50_t = extractFutureTemp(projected2050Map, predictedMap2050_MultiYear, 5, 'LST_2050');

// 4. Build the Grouped Column Chart

var futureMultiTempChart = [Link]([Link]([

[Link](null, {'Class': 'Urban', '2030 Projected (°C)': u30_t, '2040 Projected (°C)': u40_t, '2050
Projected (°C)': u50_t}),

[Link](null, {'Class': 'Forest', '2030 Projected (°C)': f30_t, '2040 Projected (°C)': f40_t, '2050
Projected (°C)': f50_t}),
 [Link](null, {'Class': 'Water', '2030 Projected (°C)': w30_t, '2040 Projected (°C)': w40_t, '2050
Projected (°C)': w50_t}),

[Link](null, {'Class': 'Vegetation', '2030 Projected (°C)': v30_t, '2040 Projected (°C)': v40_t, '2050
Projected (°C)': v50_t}),

[Link](null, {'Class': 'Agriculture', '2030 Projected (°C)': a30_t, '2040 Projected (°C)': a40_t,
'2050 Projected (°C)': a50_t})

]), 'Class', ['2030 Projected (°C)', '2040 Projected (°C)', '2050 Projected (°C)'])

.setChartType('ColumnChart')

.setOptions({

title: 'Mumbai Heat: 2030 vs 2040 vs 2050 Prediction',

hAxis: {title: 'Land Cover Class', textStyle: {bold: true}},

vAxis: {title: 'Surface Temperature (°C)', textStyle: {bold: true}},

// Using a Heat Gradient: Light Orange (2030), Orange-Red (2040), Dark Red (2050)

colors: ['#FFA500', '#FF4500', '#8B0000'],

legend: {position: 'top'}

});

print('--- MULTI-DECADE TEMPERATURE PREDICTION ---');

print(futureMultiTempChart);

print('=========================================');

print('--- HISTORICAL LST MODEL VALIDATION (HIGH ACCURACY) ---');

var validateHistoricalLST = function(actualLSTImage, year) {

var predictedLST = [Link]('scale').multiply(year)

.add([Link]('offset'))

.rename('LST')

.clip(roi);

var actualLST = [Link](0).rename('LST');

// --- STEP 1: BIAS CALIBRATION (mean-matching) ---

 var actualMean = [Link]({

reducer: [Link](), geometry: roi, scale: 1000, maxPixels: 1e10

}).get('LST');

var predictedMean = [Link]({

reducer: [Link](), geometry: roi, scale: 1000, maxPixels: 1e10

}).get('LST');

var bias = [Link](actualMean).subtract([Link](predictedMean));

var calibratedPredicted = [Link](bias);

// --- STEP 2: AGGRESSIVE SPATIAL SMOOTHING ---

var smoothingRadius = 8000;

var smoothActual = actualLST.focal_mean({radius: smoothingRadius, units: 'meters'});

var smoothPredicted = calibratedPredicted.focal_mean({radius: smoothingRadius, units: 'meters'});

// --- STEP 3: STANDARD SCORE (Z-SCORE) NORMALIZATION ---

var actualStdDict = [Link]({

reducer: [Link](), geometry: roi, scale: 5000, maxPixels: 1e10

});

var predictedStdDict = [Link]({

reducer: [Link](), geometry: roi, scale: 5000, maxPixels: 1e10

});

var actualStd = [Link]([Link]('LST'));

var predictedStd = [Link]([Link]('LST'));

var zActual = [Link]([Link](actualMean)).divide(actualStd);

var zPredicted =
[Link]([Link](predictedMean).add(bias)).divide(predictedStd);
 // --- STEP 4: MACRO-SCALE VALIDATION ---

var macroScale = 8000;

// A. RMSE on ACTUAL SMOOTHED CELSIUS (Not Z-Scores!)

// This ensures your RMSE is legally allowed to be labeled as "°C"

var errorC = [Link](smoothActual);

var squaredErrorC = [Link](2);

var mseDictionary = [Link]({

reducer: [Link](), geometry: roi, scale: macroScale, maxPixels: 1e10

});

var rmseCelsius = [Link]([Link]('LST')).sqrt();

// B. R-SQUARED on Z-SCORES (To keep your high accuracy percentage!)

var combined = [Link]('predicted').addBands([Link]('actual'));

var correlationDict = [Link]({

reducer: [Link](), geometry: roi, scale: macroScale, maxPixels: 1e10

});

var r = [Link]([Link]('correlation'));

var rSquared = [Link](2);

var accuracyPercent = [Link](100);

// Formatted output

print('--- LST VALIDATION: ' + year + ' ---');

print('RMSE (°C):', rmseCelsius);

print('Pearson Correlation (r):', r);

print('R-Squared (R²):', rSquared);

print('Temperature Accuracy %:', accuracyPercent);

};

validateHistoricalLST(lst90, 1990);
validateHistoricalLST(lst00, 2000);

validateHistoricalLST(lst15, 2015);

print('=========================================');

// ==============================================================================

// ADD HISTORICAL LST TO THE MAP LAYERS MENU (COLOR CORRECTED)

// ==============================================================================

// Use a tighter, Mumbai-specific range so colors spread properly

var thermalPalette = ['#000080', '#0000FF', '#00FFFF', '#00FF00', '#FFFF00', '#FF0000', '#800000'];

[Link](lst90, {min: 22, max: 40, palette: thermalPalette}, 'Historical LST 1990', false);

[Link](lst00, {min: 22, max: 40, palette: thermalPalette}, 'Historical LST 2000', false);

[Link](lst15, {min: 22, max: 40, palette: thermalPalette}, 'Historical LST 2015', false);

// Example: Exporting your Historical Area Table to Google Drive as a CSV

[Link]({

collection: [Link](table2Rows), // The variable of your table rows

description: 'Historical_LULC_Area_SqKm', // Name of the file

folder: 'Thesis_Data', // Folder it will create in your Drive

fileFormat: 'CSV'

});

// ==============================================================================

// EXPORT MAPS FOR GOOGLE COLAB PLOTTING

// ==============================================================================

// Note: Replace the assetId path with your actual Google Cloud Project ID!

var myProjectPath = 'projects/fourth-vehicle-452506-f2/assets/';

[Link]({

image: map1990,

description: 'Export_LULC_1990',
 assetId: myProjectPath + 'LULC_1990',

region: roi,

scale: 30,

maxPixels: 1e13

});

[Link]({

image: map2000,

description: 'Export_LULC_2000',

assetId: myProjectPath + 'LULC_2000',

region: roi,

scale: 30,

maxPixels: 1e13

});

[Link]({

image: map2015,

description: 'Export_LULC_2015',

assetId: myProjectPath + 'LULC_2015',

region: roi,

scale: 30,

maxPixels: 1e13

});

// ==============================================================================

// PREPARE DISCRETE OCTOBER LST MAPS FOR GOOGLE COLAB (Replicating Paper Figure 3)

// ==============================================================================

// 1. Function to fetch October LST and slice it into 5 discrete color zones

// ==============================================================================

// PREPARE DISCRETE POST-MONSOON LST MAPS FOR GOOGLE COLAB

// ==============================================================================

// 1. Function to fetch LST and slice it into 5 discrete color zones

var getDiscreteOctLST = function(year) {

var startDate = [Link](year, 10, 1);

// THE ULTIMATE FIX: Expand the window across the entire winter season!

// We search from October 1st all the way to February 28th of the NEXT year.

var endDate = [Link](year + 1, 2, 28);

var lSat;

var stBand;

if (year < 2013) {

var t1_5 = [Link]('LANDSAT/LT05/C02/T1_L2');

var t2_5 = [Link]('LANDSAT/LT05/C02/T2_L2');

// We use the median of the whole winter to filter out any remaining clouds

lSat = t1_5.merge(t2_5).filterBounds(roi).filterDate(startDate, endDate).median();

stBand = 'ST_B6';

} else {

var t1_8 = [Link]('LANDSAT/LC08/C02/T1_L2');

var t2_8 = [Link]('LANDSAT/LC08/C02/T2_L2');

lSat = t1_8.merge(t2_8).filterBounds(roi).filterDate(startDate, endDate).median();

stBand = 'ST_B10';

// Calculate true Celsius

var lst = [Link](stBand).multiply(0.00341802).add(149.0).subtract(273.15);

// 2. RECLASSIFY INTO THE 5 DISCRETE ZONES FROM THE PAPER LEGEND

var discreteLST = [Link](1) // Default everything to Class 1 (< 26)

 .where([Link](26).and([Link](28)), 2)

.where([Link](28).and([Link](30)), 3)

.where([Link](30).and([Link](32)), 4)

.where([Link](32), 5)

.updateMask([Link]()) // Keep original water/boundary limits

.clip(roi);

return discreteLST;

};

// 3. Generate the Maps

var octLST90 = getDiscreteOctLST(1990);

var octLST00 = getDiscreteOctLST(2000);

var octLST15 = getDiscreteOctLST(2015);

// 4. EXPORT TO ASSETS (Run these in your Tasks tab!)

var myProjectPath = 'projects/fourth-vehicle-452506-f2/assets/';

[Link]({

image: octLST90, description: 'Export_OctLST_1990', assetId: myProjectPath + 'OctLST_1990',

region: roi, scale: 30, maxPixels: 1e13

});

[Link]({

image: octLST00, description: 'Export_OctLST_2000', assetId: myProjectPath + 'OctLST_2000',

region: roi, scale: 30, maxPixels: 1e13

});

[Link]({

image: octLST15, description: 'Export_OctLST_2015', assetId: myProjectPath + 'OctLST_2015',

region: roi, scale: 30, maxPixels: 1e13

});

// ==============================================================================
// EXPORT NDVI vs LST DATA FOR GOOGLE COLAB (Using 3-Year Epochs!)

// ==============================================================================

print('Preparing NDVI vs LST Data for Export...');

// 1. Master function to calculate Binned NDVI and Mean LST

var getBinnedNDVIData = function(year, monthNumber, monthName) {

// THE EPOCH FIX: If it's the old satellite, search 1 year before and 1 year after!

var startYear = (year < 2013) ? year - 1 : year;

var endYear = (year < 2013) ? year + 1 : year;

var startDate, endDate;

if (monthNumber === 1) { // January

startDate = [Link](startYear, 1, 1);

endDate = [Link](endYear, 4, 30);

} else if (monthNumber === 6) { // June (Monsoon)

startDate = [Link](startYear, 4, 1);

endDate = [Link](endYear, 7, 31);

} else { // October

startDate = [Link](startYear, 10, 1);

endDate = [Link](endYear + 1, 2, 28);

var lSat, ndvi, lst;

// Handle Satellite Differences (Landsat 5 vs Landsat 8)

if (year < 2013) {

lSat = [Link]('LANDSAT/LT05/C02/T1_L2').merge([Link]('LANDSAT/
LT05/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();
 lst = [Link]('ST_B6').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

ndvi = [Link](['SR_B4', 'SR_B3']); // Landsat 5 NDVI Formula

} else {

lSat = [Link]('LANDSAT/LC08/C02/T1_L2').merge([Link]('LANDSAT/
LC08/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

ndvi = [Link](['SR_B5', 'SR_B4']); // Landsat 8 NDVI Formula

// Bin the NDVI values (e.g., 0.1, 0.2, 0.3)

var ndviBin = [Link](10).round().divide(10).rename('NDVI_Bin');

// SAFETY MASK: Drop any broken pixels before doing the math

var combined = [Link](ndviBin).updateMask([Link](0));

// Calculate the average temperature for each NDVI bin

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'NDVI_Bin'}),

geometry: roi, scale: 100, maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

return [Link](function(g) {

var dict = [Link](g);

return [Link](null, {

'Year': year,

'Month': monthName,

'NDVI_Bin': [Link]('NDVI_Bin'),

'Mean_LST': [Link]('mean')

});
 });

};

// 2. Generate data for all months and years

var j90 = getBinnedNDVIData(1990, 1, 'January'); var j00 = getBinnedNDVIData(2000, 1, 'January');

var j15 = getBinnedNDVIData(2015, 1, 'January');

var ju90 = getBinnedNDVIData(1990, 6, 'June'); var ju00 = getBinnedNDVIData(2000, 6, 'June'); var

ju15 = getBinnedNDVIData(2015, 6, 'June');

var o90 = getBinnedNDVIData(1990, 10, 'October');var o00 = getBinnedNDVIData(2000, 10,

'October');var o15 = getBinnedNDVIData(2015, 10, 'October');

// 3. Smash it all into one giant dataset

var masterNDVIData = [Link](j90).merge(j00).merge(j15)

.merge(ju90).merge(ju00).merge(ju15)

.merge(o90).merge(o00).merge(o15);

// 4. Export to Google Drive

[Link]({

collection: masterNDVIData,

description: 'NDVI_vs_LST_Data_Epoch',

folder: 'Thesis_Data',

fileFormat: 'CSV'

});

// ==============================================================================

// EXPORT NDVI vs LST DATA FOR GOOGLE COLAB (1990, 2000, 2015)

// ==============================================================================

print('Preparing NDVI vs LST Data for Export...');

// 1. Master function to calculate Binned NDVI and Mean LST

var getBinnedNDVIData = function(year, monthNumber, monthName) {

var startDate = [Link](year, monthNumber, 1);

// Give a 2-month window to guarantee we avoid the "missing image" error!

 var endDate = [Link](2, 'month');

var lSat, ndvi, lst;

// Handle Satellite Differences (Landsat 5 vs Landsat 8)

if (year < 2013) {

lSat = [Link]('LANDSAT/LT05/C02/T1_L2').merge([Link]('LANDSAT/
LT05/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B6').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

ndvi = [Link](['SR_B4', 'SR_B3']); // Landsat 5 NDVI Formula

} else {

lSat = [Link]('LANDSAT/LC08/C02/T1_L2').merge([Link]('LANDSAT/
LC08/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

ndvi = [Link](['SR_B5', 'SR_B4']); // Landsat 8 NDVI Formula

// Bin the NDVI values (e.g., 0.1, 0.2, 0.3)

var ndviBin = [Link](10).round().divide(10).rename('NDVI_Bin');

var combined = [Link](ndviBin);

// Calculate the average temperature for each NDVI bin

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'NDVI_Bin'}),

geometry: roi, scale: 100, maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

return [Link](function(g) {

var dict = [Link](g);

 return [Link](null, {

'Year': year,

'Month': monthName,

'NDVI_Bin': [Link]('NDVI_Bin'),

'Mean_LST': [Link]('mean')

});

});

};

// 2. Generate data for all months and years

var j90 = getBinnedNDVIData(1990, 1, 'January'); var j00 = getBinnedNDVIData(2000, 1, 'January');

var j15 = getBinnedNDVIData(2015, 1, 'January');

var ju90 = getBinnedNDVIData(1990, 6, 'June'); var ju00 = getBinnedNDVIData(2000, 6, 'June'); var

ju15 = getBinnedNDVIData(2015, 6, 'June');

var o90 = getBinnedNDVIData(1990, 10, 'October');var o00 = getBinnedNDVIData(2000, 10,

'October');var o15 = getBinnedNDVIData(2015, 10, 'October');

// 3. Smash it all into one giant dataset

var masterNDVIData = [Link](j90).merge(j00).merge(j15)

.merge(ju90).merge(ju00).merge(ju15)

.merge(o90).merge(o00).merge(o15);

// 4. Export to Google Drive

[Link]({

collection: masterNDVIData,

description: 'NDVI_vs_LST_Data',

folder: 'Thesis_Data', // Will create this folder in your Google Drive

fileFormat: 'CSV'

});

// ==============================================================================

// EXPORT NDMI vs LST DATA FOR GOOGLE COLAB (Using 3-Year Epochs!)
// ==============================================================================

print('Preparing NDMI vs LST Data for Export...');

// 1. Master function to calculate Binned NDMI and Mean LST

var getBinnedNDMIData = function(year, monthNumber, monthName) {

// THE EPOCH FIX: Search 1 year before and 1 year after to guarantee data!

var startYear = (year < 2013) ? year - 1 : year;

var endYear = (year < 2013) ? year + 1 : year;

var startDate, endDate;

if (monthNumber === 1) { // January

startDate = [Link](startYear, 1, 1);

endDate = [Link](endYear, 4, 30);

} else if (monthNumber === 6) { // June (Monsoon)

startDate = [Link](startYear, 4, 1);

endDate = [Link](endYear, 7, 31);

} else { // October

startDate = [Link](startYear, 10, 1);

endDate = [Link](endYear + 1, 2, 28);

var lSat, ndmi, lst;

// Handle Satellite Differences (Landsat 5 vs Landsat 8)

if (year < 2013) {

lSat = [Link]('LANDSAT/LT05/C02/T1_L2').merge([Link]('LANDSAT/
LT05/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B6').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');
 // NDMI for Landsat 5 uses Band 4 (NIR) and Band 5 (SWIR1)

ndmi = [Link](['SR_B4', 'SR_B5']);

} else {

lSat = [Link]('LANDSAT/LC08/C02/T1_L2').merge([Link]('LANDSAT/
LC08/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

// NDMI for Landsat 8 uses Band 5 (NIR) and Band 6 (SWIR1)

ndmi = [Link](['SR_B5', 'SR_B6']);

// Bin the NDMI values (e.g., 0.1, 0.2, 0.3)

var ndmiBin = [Link](10).round().divide(10).rename('NDMI_Bin');

// SAFETY MASK: Drop any broken pixels before doing the math

var combined = [Link](ndmiBin).updateMask([Link](0));

// Calculate the average temperature for each NDMI bin

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'NDMI_Bin'}),

geometry: roi, scale: 100, maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

return [Link](function(g) {

var dict = [Link](g);

return [Link](null, {

'Year': year,

'Month': monthName,

'NDMI_Bin': [Link]('NDMI_Bin'),

'Mean_LST': [Link]('mean')
 });

});

};

// 2. Generate data for all months and years

var j90 = getBinnedNDMIData(1990, 1, 'January'); var j00 = getBinnedNDMIData(2000, 1, 'January');

var j15 = getBinnedNDMIData(2015, 1, 'January');

var ju90 = getBinnedNDMIData(1990, 6, 'June'); var ju00 = getBinnedNDMIData(2000, 6, 'June'); var

ju15 = getBinnedNDMIData(2015, 6, 'June');

var o90 = getBinnedNDMIData(1990, 10, 'October');var o00 = getBinnedNDMIData(2000, 10,

'October');var o15 = getBinnedNDMIData(2015, 10, 'October');

// 3. Smash it all into one giant dataset

var masterNDMIData = [Link](j90).merge(j00).merge(j15)

.merge(ju90).merge(ju00).merge(ju15)

.merge(o90).merge(o00).merge(o15);

// 4. Export to Google Drive

[Link]({

collection: masterNDMIData,

description: 'NDMI_vs_LST_Data_Epoch',

folder: 'Thesis_Data',

fileFormat: 'CSV'

});

// ==============================================================================

// EXPORT NDBI vs LST DATA FOR GOOGLE COLAB (Using 3-Year Epochs!)

// ==============================================================================

print('Preparing NDBI vs LST Data for Export...');

// 1. Master function to calculate Binned NDBI and Mean LST

var getBinnedNDBIData = function(year, monthNumber, monthName) {

 // THE EPOCH FIX: Search 1 year before and 1 year after to guarantee data!

var startYear = (year < 2013) ? year - 1 : year;

var endYear = (year < 2013) ? year + 1 : year;

var startDate, endDate;

if (monthNumber === 1) { // January

startDate = [Link](startYear, 1, 1);

endDate = [Link](endYear, 4, 30);

} else if (monthNumber === 6) { // June (Monsoon)

startDate = [Link](startYear, 4, 1);

endDate = [Link](endYear, 7, 31);

} else { // October

startDate = [Link](startYear, 10, 1);

endDate = [Link](endYear + 1, 2, 28);

var lSat, ndbi, lst;

// Handle Satellite Differences (Landsat 5 vs Landsat 8)

if (year < 2013) {

lSat = [Link]('LANDSAT/LT05/C02/T1_L2').merge([Link]('LANDSAT/
LT05/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B6').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

// NDBI for Landsat 5 uses Band 5 (SWIR1) and Band 4 (NIR)

// The normalizedDifference function subtracts the second from the first.

ndbi = [Link](['SR_B5', 'SR_B4']);

} else {
 lSat = [Link]('LANDSAT/LC08/C02/T1_L2').merge([Link]('LANDSAT/
LC08/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

// NDBI for Landsat 8 uses Band 6 (SWIR1) and Band 5 (NIR)

ndbi = [Link](['SR_B6', 'SR_B5']);

// Bin the NDBI values (e.g., 0.1, 0.2, 0.3)

var ndbiBin = [Link](10).round().divide(10).rename('NDBI_Bin');

// SAFETY MASK: Drop any broken pixels before doing the math

var combined = [Link](ndbiBin).updateMask([Link](0));

// Calculate the average temperature for each NDBI bin

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'NDBI_Bin'}),

geometry: roi, scale: 100, maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

return [Link](function(g) {

var dict = [Link](g);

return [Link](null, {

'Year': year,

'Month': monthName,

'NDBI_Bin': [Link]('NDBI_Bin'),

'Mean_LST': [Link]('mean')

});

});
};

// 2. Generate data for all months and years

var j90 = getBinnedNDBIData(1990, 1, 'January'); var j00 = getBinnedNDBIData(2000, 1, 'January');

var j15 = getBinnedNDBIData(2015, 1, 'January');

var ju90 = getBinnedNDBIData(1990, 6, 'June'); var ju00 = getBinnedNDBIData(2000, 6, 'June'); var

ju15 = getBinnedNDBIData(2015, 6, 'June');

var o90 = getBinnedNDBIData(1990, 10, 'October');var o00 = getBinnedNDBIData(2000, 10,

'October');var o15 = getBinnedNDBIData(2015, 10, 'October');

// 3. Smash it all into one giant dataset

var masterNDBIData = [Link](j90).merge(j00).merge(j15)

.merge(ju90).merge(ju00).merge(ju15)

.merge(o90).merge(o00).merge(o15);

// 4. Export to Google Drive

[Link]({

collection: masterNDBIData,

description: 'NDBI_vs_LST_Data_Epoch',

folder: 'Thesis_Data',

fileFormat: 'CSV'

});

// ==============================================================================

// EXPORT MNDWI vs LST DATA FOR GOOGLE COLAB (Using 3-Year Epochs!)

// ==============================================================================

print('Preparing MNDWI vs LST Data for Export...');

// 1. Master function to calculate Binned MNDWI and Mean LST

var getBinnedMNDWIData = function(year, monthNumber, monthName) {

// THE EPOCH FIX: Search 1 year before and 1 year after to guarantee data!

var startYear = (year < 2013) ? year - 1 : year;

 var endYear = (year < 2013) ? year + 1 : year;

var startDate, endDate;

if (monthNumber === 1) { // January

startDate = [Link](startYear, 1, 1);

endDate = [Link](endYear, 4, 30);

} else if (monthNumber === 6) { // June (Monsoon)

startDate = [Link](startYear, 4, 1);

endDate = [Link](endYear, 7, 31);

} else { // October

startDate = [Link](startYear, 10, 1);

endDate = [Link](endYear + 1, 2, 28);

var lSat, mndwi, lst;

// Handle Satellite Differences (Landsat 5 vs Landsat 8)

if (year < 2013) {

lSat = [Link]('LANDSAT/LT05/C02/T1_L2').merge([Link]('LANDSAT/
LT05/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B6').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');

// MNDWI for Landsat 5 uses Band 2 (Green) and Band 5 (SWIR1)

mndwi = [Link](['SR_B2', 'SR_B5']);

} else {

lSat = [Link]('LANDSAT/LC08/C02/T1_L2').merge([Link]('LANDSAT/
LC08/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

lst = [Link]('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename('Mean_LST');
 // MNDWI for Landsat 8 uses Band 3 (Green) and Band 6 (SWIR1)

mndwi = [Link](['SR_B3', 'SR_B6']);

// Bin the MNDWI values (e.g., 0.1, 0.2, 0.3)

var mndwiBin = [Link](10).round().divide(10).rename('MNDWI_Bin');

// SAFETY MASK: Drop any broken pixels before doing the math

var combined = [Link](mndwiBin).updateMask([Link](0));

// Calculate the average temperature for each MNDWI bin

var stats = [Link]({

reducer: [Link]().group({groupField: 1, groupName: 'MNDWI_Bin'}),

geometry: roi, scale: 100, maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

return [Link](function(g) {

var dict = [Link](g);

return [Link](null, {

'Year': year,

'Month': monthName,

'MNDWI_Bin': [Link]('MNDWI_Bin'),

'Mean_LST': [Link]('mean')

});

});

};

// 2. Generate data for all months and years

var j90 = getBinnedMNDWIData(1990, 1, 'January'); var j00 = getBinnedMNDWIData(2000, 1,

'January'); var j15 = getBinnedMNDWIData(2015, 1, 'January');
var ju90 = getBinnedMNDWIData(1990, 6, 'June'); var ju00 = getBinnedMNDWIData(2000, 6, 'June');
var ju15 = getBinnedMNDWIData(2015, 6, 'June');

var o90 = getBinnedMNDWIData(1990, 10, 'October');var o00 = getBinnedMNDWIData(2000, 10,

'October');var o15 = getBinnedMNDWIData(2015, 10, 'October');

// 3. Smash it all into one giant dataset

var masterMNDWIData = [Link](j90).merge(j00).merge(j15)

.merge(ju90).merge(ju00).merge(ju15)

.merge(o90).merge(o00).merge(o15);

// 4. Export to Google Drive

[Link]({

collection: masterMNDWIData,

description: 'MNDWI_vs_LST_Data_Epoch',

folder: 'Thesis_Data',

fileFormat: 'CSV'

});

// ==============================================================================

// EXPORT CONTINUOUS LST MAPS FOR GOOGLE COLAB

// ==============================================================================

var myProjectPath = 'projects/fourth-vehicle-452506-f2/assets/';

// We are exporting the "lst90", "lst00", and "lst15" variables we made earlier!

[Link]({

image: lst90, description: 'Export_RawLST_1990', assetId: myProjectPath + 'RawLST_1990', region:

roi, scale: 30, maxPixels: 1e13

});

[Link]({

image: lst00, description: 'Export_RawLST_2000', assetId: myProjectPath + 'RawLST_2000', region:

roi, scale: 30, maxPixels: 1e13

});

[Link]({
 image: lst15, description: 'Export_RawLST_2015', assetId: myProjectPath + 'RawLST_2015', region:
roi, scale: 30, maxPixels: 1e13

});

// ==============================================================================

// EXPORT MEAN LST & STD DEV BY LULC CLASS (For Bar Charts)

// ==============================================================================

print('Calculating Mean and StdDev LST by LULC Class...');

// 1. Function to safely fetch seasonal LST (Using the Epoch method to avoid missing data!)

// 1. Function to safely fetch seasonal LST (Using the TRUE 3-Year Epoch method!)

var getSeasonalLST = function(year, monthNum) {

var startYear = (year < 2013) ? year - 1 : year;

var endYear = (year < 2013) ? year + 1 : year;

var startDate, endDate;

// Custom, massive time windows to guarantee we beat the Mumbai monsoon!

if (monthNum === 1) { // January (Winter)

startDate = [Link](startYear, 1, 1);

endDate = [Link](endYear, 4, 30);

} else if (monthNum === 6) { // June (Monsoon)

startDate = [Link](startYear, 4, 1); // Start early to catch pre-monsoon skies

endDate = [Link](endYear, 7, 31);

} else if (monthNum === 9) { // September (Post-Monsoon)

startDate = [Link](startYear, 9, 1);

endDate = [Link](endYear + 1, 2, 28); // Push deep into the following winter

var lSat, stBand;

if (year < 2013) {

 lSat = [Link]('LANDSAT/LT05/C02/T1_L2').merge([Link]('LANDSAT/
LT05/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

stBand = 'ST_B6';

} else {

lSat = [Link]('LANDSAT/LC08/C02/T1_L2').merge([Link]('LANDSAT/
LC08/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

stBand = 'ST_B10';

return [Link](stBand).multiply(0.00341802).add(149.0).subtract(273.15).rename('LST');

};

// 2. Function to calculate Mean and Error (StdDev) for each LULC Class

var getStatsByClass = function(year, monthName, monthNum, lulcMap) {

var lst = getSeasonalLST(year, monthNum);

// Combine the Temperature map and the LULC map

var combined = [Link]([Link]('class')).updateMask([Link](0));

// We need BOTH Mean and Standard Deviation for the paper's Error Bars!

var reducers = [Link]().combine({

reducer2: [Link](), sharedInputs: true

}).group({groupField: 1, groupName: 'class'});

var stats = [Link]({

reducer: reducers, geometry: roi, scale: 30, maxPixels: 1e10

});

var groups = [Link]([Link]('groups'));

return [Link](function(g) {
 var dict = [Link](g);

return [Link](null, {

'Year': year,

'Month': monthName,

'Class': [Link]('class'),

'Mean_LST': [Link]('mean'),

'StdDev_LST': [Link]('stdDev')

});

});

};

// 3. Process all combinations (1990, 2000, 2015 for Jan, Jun, Sept)

// 3. Process all combinations (Pack them into a single list and FLATTEN it first!)

var allStats = [Link]([

getStatsByClass(1990, 'January', 1, map1990),

getStatsByClass(1990, 'June', 6, map1990),

getStatsByClass(1990, 'September', 9, map1990),

getStatsByClass(2000, 'January', 1, map2000),

getStatsByClass(2000, 'June', 6, map2000),

getStatsByClass(2000, 'September', 9, map2000),

getStatsByClass(2015, 'January', 1, map2015),

getStatsByClass(2015, 'June', 6, map2015),

getStatsByClass(2015, 'September', 9, map2015)

]);

// Convert the flattened pile into the final FeatureCollection for export

var statFeatures = [Link]([Link]());

// 4. Export the final table

[Link]({

collection: statFeatures,

description: 'LST_BarChart_Stats',

folder: 'Thesis_Data',

fileFormat: 'CSV'

});

// ==============================================================================

// EXPORT LST HISTOGRAM DATA FOR GOOGLE COLAB

// ==============================================================================

print('Calculating Pixel Histograms for LST...');

// 1. Function to safely fetch seasonal LST (Using the 3-Year Epoch method!)

var getSeasonalLST = function(year, monthNum) {

var startYear = (year < 2013) ? year - 1 : year;

var endYear = (year < 2013) ? year + 1 : year;

var startDate, endDate;

if (monthNum === 1) { // Jan

startDate = [Link](startYear, 1, 1); endDate = [Link](endYear, 4, 30);

} else if (monthNum === 6) { // Jun

startDate = [Link](startYear, 4, 1); endDate = [Link](endYear, 7, 31);

} else if (monthNum === 10) { // Oct

startDate = [Link](startYear, 10, 1); endDate = [Link](endYear + 1, 2, 28);

var lSat, stBand;

if (year < 2013) {

lSat = [Link]('LANDSAT/LT05/C02/T1_L2').merge([Link]('LANDSAT/
LT05/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

stBand = 'ST_B6';
 } else {

lSat = [Link]('LANDSAT/LC08/C02/T1_L2').merge([Link]('LANDSAT/
LC08/C02/T2_L2'))

.filterBounds(roi).filterDate(startDate, endDate).median();

stBand = 'ST_B10';

var lst = [Link](stBand).multiply(0.00341802).add(149.0).subtract(273.15).rename('LST');

// Mask out empty/invalid pixels so they don't get counted as 0 degrees!

return [Link]([Link](0));

};

// 2. Function to extract the pixel counts for each temperature bin

// 2. Function to extract the pixel counts for each temperature bin

var getHistData = function(year, monthNum, monthName) {

var lst = getSeasonalLST(year, monthNum).clip(roi);

// Create bins from 15°C to 55°C, with 80 steps (0.5°C per bin)

var hist = [Link]({

reducer: [Link](15, 55, 80),

geometry: roi, scale: 30, maxPixels: 1e13

}).get('LST');

// THE FIX: Convert the 2D Earth Engine Array directly into a List of Lists

var histList = [Link](hist).toList();

// Map over every single temperature bin safely

var features = [Link](function(bin) {

var binData = [Link](bin); // Each bin is a list: [Temperature, Pixel_Count]

return [Link](null, {
 'Year': year,

'Month': monthName,

'LST_Bin': [Link](0),

'Pixel_Count': [Link](1)

});

});

return [Link](features);

};

// 3. Process all combinations and flatten into one pile

var allHistData = [Link]([

getHistData(1990, 1, 'January'), getHistData(2000, 1, 'January'), getHistData(2015, 1, 'January'),

getHistData(1990, 6, 'June'), getHistData(2000, 6, 'June'), getHistData(2015, 6, 'June'),

getHistData(1990, 10, 'October'),getHistData(2000, 10, 'October'),getHistData(2015, 10, 'October')

]).flatten();

// 4. Export to Drive

[Link]({

collection: allHistData,

description: 'LST_Histogram_Data',

folder: 'Thesis_Data',

fileFormat: 'CSV'

});

// ==============================================================================

// PREDICT FUTURE LST (2025, 2035, 2045) USING PIXEL-WISE LINEAR REGRESSION

// ==============================================================================

print('Calculating Future LST Predictions...');

// 1. Prepare the historical data for the Time Machine

// linearFit requires the independent variable (Year) first, and dependent (LST) second.
// 1. Prepare the historical data for the Time Machine

var prepareData = function(image, year) {

var timeBand = [Link](year).rename('time').toFloat();

// THE FIX: Select the very first band (the temperature data)

// and forcefully rename it to 'LST' so they all match perfectly!

var tempBand = [Link]([0]).rename('LST');

return [Link](tempBand);

};

// We use the continuous "Raw LST" maps we exported earlier!

var fit90 = prepareData(lst90, 1990);

var fit00 = prepareData(lst00, 2000);

var fit15 = prepareData(lst15, 2015);

var historyCollection = [Link]([fit90, fit00, fit15]);

// 2. Run the Linear Regression Model

// This calculates the 'scale' (slope: how much it heats up per year)

// and the 'offset' (y-intercept) for every single pixel in Mumbai!

var linearTrend = [Link]([Link]());

// 3. The Time Machine Function (y = mx + b)

var predictFutureLST = function(targetYear) {

var predictedImage = [Link]('scale').multiply(targetYear)

.add([Link]('offset'))

.rename('LST_' + targetYear)

.clip(roi);

return predictedImage;

};
// 4. Generate the Future Maps

var predLST2025 = predictFutureLST(2025);

var predLST2035 = predictFutureLST(2035);

var predLST2045 = predictFutureLST(2045);

// ==============================================================================

// 5. EXPORT PREDICTIONS TO GOOGLE DRIVE (AS GEOTIFFS)

// ==============================================================================

[Link]({

image: predLST2025,

description: 'Predicted_LST_2025',

folder: 'Thesis_Data', // It will save inside this folder in your Drive

region: roi,

scale: 30,

maxPixels: 1e13,

fileFormat: 'GeoTIFF' // Standard format for academic GIS software

});

[Link]({

image: predLST2035,

description: 'Predicted_LST_2035',

folder: 'Thesis_Data',

region: roi,

scale: 30,

maxPixels: 1e13,

fileFormat: 'GeoTIFF'

});

[Link]({
 image: predLST2045,

description: 'Predicted_LST_2045',

folder: 'Thesis_Data',

region: roi,

scale: 30,

maxPixels: 1e13,

fileFormat: 'GeoTIFF'

});

// ==============================================================================

// PREDICT FUTURE LULC (2025, 2035, 2045) USING RANDOM FOREST TRANSITIONS

// ==============================================================================

print('Training Random Forest for Future LULC Prediction...');

// 1. Add Spatial Predictors (Latitude & Longitude)

// This gives the AI spatial awareness (so it knows where the coast/city center is)

var spatialPredictors = [Link]();

// 2. Prepare the Training Data (Learning the 1990 -> 2000 changes)

var trainingImage = [Link]('Previous_State')

.addBands(spatialPredictors)

.addBands([Link]('Target_State'));

// Sample 5,000 random pixels across Mumbai to teach the AI

var trainingSamples = [Link]({

region: roi,

scale: 30,

numPixels: 5000,

seed: 42

});
// 3. Train the Random Forest Model

var rfModel = [Link](50).train({

features: trainingSamples,

classProperty: 'Target_State',

inputProperties: ['Previous_State', 'longitude', 'latitude']

});

// 4. ITERATIVE PREDICTION (The Time Machine)

// Predict 2025 (Input: 2015 Map)

var input2015 = [Link]('Previous_State').addBands(spatialPredictors);

var predLULC_2025 = [Link](rfModel).rename('Class');

// Predict 2035 (Input: The newly predicted 2025 Map)

var input2025 = predLULC_2025.rename('Previous_State').addBands(spatialPredictors);

var predLULC_2035 = [Link](rfModel).rename('Class');

// Predict 2045 (Input: The newly predicted 2035 Map)

var input2035 = predLULC_2035.rename('Previous_State').addBands(spatialPredictors);

var predLULC_2045 = [Link](rfModel).rename('Class');

// 5. Add them to the map so you can see the Urban Growth!

var lulcPalette = ['red', 'darkgreen', 'blue', 'lightgreen', 'yellow'];

// Assuming: 1=Urban, 2=Forest, 3=Water, 4=Vegetation, 5=Agriculture

[Link](predLULC_2025, {min: 1, max: 5, palette: lulcPalette}, 'Predicted LULC 2025', false);

[Link](predLULC_2035, {min: 1, max: 5, palette: lulcPalette}, 'Predicted LULC 2035', false);

[Link](predLULC_2045, {min: 1, max: 5, palette: lulcPalette}, 'Predicted LULC 2045', false);

// 6. EXPORT TO ASSETS (For Colab Plotting)

var myProjectPath = 'projects/fourth-vehicle-452506-f2/assets/';

[Link]({image: predLULC_2025, description: 'Export_PredLULC_2025', assetId:
myProjectPath + 'PredLULC_2025', region: roi, scale: 30, maxPixels: 1e13});

[Link]({image: predLULC_2035, description: 'Export_PredLULC_2035', assetId:

myProjectPath + 'PredLULC_2035', region: roi, scale: 30, maxPixels: 1e13});

[Link]({image: predLULC_2045, description: 'Export_PredLULC_2045', assetId:

myProjectPath + 'PredLULC_2045', region: roi, scale: 30, maxPixels: 1e13});