// A wall–to–wall comparison of MapBiomas (Level 1) and JRC–Tropical Moist Forests 
// Calculating exact transition areas in square meters (m2) across the entire MapBiomas Amazon footprint
   
//######################################################################
// SECTION 1 : IMPORT BASE DATASETS & PALETTES (2023)
//######################################################################


// Load JRC TMF Annual Changes Collection & Mosaic : Select year 2024
var AnnualChanges = ee.ImageCollection('projects/JRC/TMF/v1_2025/AnnualChanges').mosaic();
var TMF = AnnualChanges.select('Dec2023');


// Load MapBiomas - Collection 6 (Pan-Amazon Asset)
var mapbiomasAmazonImage = ee.Image('projects/mapbiomas-public/assets/amazon/lulc/collection6/mapbiomas_collection60_integration_v1');


// Select year 2023 classification
var mapbiomas2023 = mapbiomasAmazonImage.select('classification_2023');


// Extract the exact ecological mapping footprint of MapBiomas
var amazonEcologicalRegion = mapbiomas2023.geometry();


//######################################################################
// SECTION 2 : RECLASSIFICATION TO LEVEL 1, ALIGNMENT & AREA SUMMATION
//######################################################################


// 1. Isolate single bands
var s1_tmf_raw = TMF.select([0]).rename('JRC_TMF_Class');
var s2_mb_raw  = mapbiomas2023.select([0]); 


// 2. RECLASSIFICATION: Mapping Level 2 Subclasses to official Level 1 IDs
var fromValues = [
  3, 4, 5, 6,             // Subclasses for Main Class 1 (Natural Forest)
  11, 12, 29, 13,         // Subclasses for Main Class 2 (Non-forest natural formation)
  15, 18, 9, 35, 21,      // Subclasses for Main Class 3 (Farming and silviculture)
  23, 24, 30, 68, 25,     // Subclasses for Main Class 4 (Non-vegetated area)
  33, 34,                 // Subclasses for Main Class 26 (Water)
  27                      // Main Class 27 (Not observed)
];


var toValues = [
  1, 1, 1, 1,             // Map to 1
  2, 2, 2, 2,             // Map to 2
  3, 3, 3, 3, 3,          // Map to 3
  4, 4, 4, 4, 4,          // Map to 4
  26, 26,                 // Map to 26
  27                      // Map to 27
];


// Execute the remapping
var s2_mb_level1 = s2_mb_raw.remap(fromValues, toValues, 27).rename('MapBiomas_Level1');


// 3. PERFECT ALIGNMENT ANCHOR
var mbProj = s2_mb_level1.projection();
var s1_tmf = s1_tmf_raw.reproject({
  crs: mbProj,
  scale: mbProj.nominalScale()
}).updateMask(s2_mb_level1.mask()); 


// 4. Generate the combination code band
var combinedCodeImage = s1_tmf.multiply(1000).add(s2_mb_level1).rename('Combination_Code');


// 5. AREA GENERATION: Create an image representing the square meters of each pixel,
// then add our transition codes as a secondary band.
var areaImage = ee.Image.pixelArea().rename('area').addBands(combinedCodeImage);


// 6. Run a grouped sum reduction over the footprint
var areaStats = areaImage.reduceRegion({
  reducer: ee.Reducer.sum().group({
    groupField: 1,            // Group based on the 'Combination_Code' band (index 1)
    groupName: 'Combination_Code'
  }),
  geometry: amazonEcologicalRegion, 
  scale: 30,                     
  maxPixels: 1e13
});


//######################################################################
// SECTION 3 : CUSTOM VISUALIZATION PARAMETERS 
//######################################################################


// Visualization parameters for TMF 
function rgb(r, g, b) {  
    var bin = r << 16 | g << 8 | b;  
    return (function(h) {  
        return new Array(7 - h.length).join("0") + h;  
    })(bin.toString(16).toUpperCase());  
}


var PALETTEAnnualChanges = [
    rgb(0,90,0), rgb(100,155,35), rgb(255,135,15), rgb(210,250,60), rgb(0,140,190), rgb(255,255,255)
];


Map.addLayer(TMF.clip(amazonEcologicalRegion), {  
    'min': 1, 'max': 6, 'palette': PALETTEAnnualChanges  
}, "JRC-Annual Changes(Year:2024)", false);


// Define MapBiomas Level 1 discrete color mappings
var mbLevel1Colors = [
  '#1f8d49', // 1: Natural Forest
  '#d6bc74', // 2: Non-forest natural formation
  '#ffefc3', // 3: Farming and Sylviculture
  '#d4271e'  // 4: Non-vegetated area
];


// Since class 26 and 27 are far apart from 4, we build a specific re-mapped 
// visualization image to align perfectly with a continuous 0-5 palette index.
var mbVisImage = s2_mb_level1.remap(
  [1, 2, 3, 4, 26, 27], // Original reclassified IDs
  [0, 1, 2, 3,  4,  5], // Sequential palette index position
  5                     // Default fallback
);


var mbPaletteSeq = [
  '#1f8d49', // Index 0 -> Class 1 (Natural Forest)
  '#d6bc74', // Index 1 -> Class 2 (Non-forest)
  '#ffefc3', // Index 2 -> Class 3 (Farming)
  '#d4271e', // Index 3 -> Class 4 (Non-veg)
  '#2532e4', // Index 4 -> Class 26 (Water)
  '#f0f0f0'  // Index 5 -> Class 27 (Not observed)
];


// Add MapBiomas Level 1 to interactive Map Canvas
Map.addLayer(mbVisImage.clip(amazonEcologicalRegion), {
  'min': 0, 
  'max': 5, 
  'palette': mbPaletteSeq
}, "MapBiomas - Level 1 Macro Classes (2023)", true);




//######################################################################
// SECTION 4 : DECODE CODES AND PREPARE FOR EXPORT
//######################################################################


// Extract the grouped list of dictionaries containing codes and summed areas
var groupsList = ee.List(areaStats.get('groups'));


var finalSankeyTable = ee.FeatureCollection(groupsList.map(function(item) {
  var dict = ee.Dictionary(item);
  var code = ee.Number(dict.get('Combination_Code'));
  var areaSquareMeters = ee.Number(dict.get('sum')); // Keeps raw square meters
  
  // Unpack the 1000-multiplied key
  var tmfSource = code.divide(1000).floor();
  var mbTarget = code.mod(1000);
  
  return ee.Feature(null, {
    'Source_JRC_TMF': tmfSource,
    'Target_MapBiomas_L1': mbTarget,
    'Area_m2': areaSquareMeters
  });
}));


/*
Export.table.toDrive({
  collection: finalSankeyTable,
  description: 'TMF_MapBiomas_Level1_Area_m2',
  folder: 'Article2_Exports', 
  fileFormat: 'CSV',
  selectors: ['Source_JRC_TMF', 'Target_MapBiomas_L1', 'Area_m2']
});
*/


//######################################################################
// SECTION 5 : ILLUSTRATE 2 CASES OF DISAGREEMENT BETWEEN PRODUCTS
//######################################################################


// CREATE A LIGHTWEIGHT VISUAL ALIGNMENT IMAGE (No hard .reproject() calculation)
var s1_tmf_vis = s1_tmf_raw.updateMask(s2_mb_level1.mask());


// CASE 1: MapBiomas Forest (1) vs TMF Deforested (3) Using the Visual Tracker
var mask1_forest_vs_defor_vis = s2_mb_level1.eq(1).and(s1_tmf_vis.eq(3));


// Renders instantly without hitting the 2^31 pixel memory ceiling
Map.addLayer(mask1_forest_vs_defor_vis.selfMask(), {
  palette: ['#FFDBBB'] 
}, 'NatForest(MapB)_vs_DeforestedLand(TMF)');


// CREATE A LIGHTWEIGHT VISUAL ALIGNMENT IMAGE (No hard .reproject() calculation)
var s1_tmf_vis = s1_tmf_raw.updateMask(s2_mb_level1.mask());


// CASE 2: MapBiomas Forest (1) vs TMF Other land cover (3) Using the Visual Tracker
var mask2_forest_vs_other_vis = s2_mb_level1.eq(1).and(s1_tmf_vis.eq(6));


// Renders instantly without hitting the 2^31 pixel memory ceiling
Map.addLayer(mask2_forest_vs_other_vis.selfMask(), {
  palette: ['#90EE90'] 
}, 'NatForest(MapB)_vs_Other(TMF)');
//Map.centerObject(Xingu, 18);


//######################################################################
// SECTION 6 : TARGETED TIF EXPORTS WITH PARSABLE CRS (XINGU TILE)
//######################################################################


/*
var xinguROI = Xingu; 
// Define a standard, universally parsable CRS string for the Amazon
// 'EPSG:4326' works perfectly for unprojected lat/lon grids
var exportCRS = 'EPSG:4326'; 


// EXPORT 1: MapBiomas Level 1 Classification
Export.image.toDrive({
  image: s2_mb_level1.uint8(), 
  description: 'Xingu_MapBiomas_Level1_30m',
  folder: 'Article2_Xingu_Panels',
  region: xinguROI,
  scale: 30, 
  crs: exportCRS, // FIXED: Replaced mbProj object with explicit clean string
  maxPixels: 1e9
});


// EXPORT 2: JRC-TMF Targeted Classes
var tmfTargetedClasses = s1_tmf_raw.updateMask(s1_tmf_raw.eq(3).or(s1_tmf_raw.eq(6)));


Export.image.toDrive({
  image: tmfTargetedClasses.unmask(0).uint8(), 
  description: 'Xingu_TMF_Defor_And_Other_30m',
  folder: 'Article2_Xingu_Panels',
  region: xinguROI,
  scale: 30,
  crs: exportCRS, // FIXED: Replaced mbProj object with explicit clean string
  maxPixels: 1e9
});


// EXPORT 3: Disagreement Masks
var combinedDisagreementMasks = mask1_forest_vs_defor_vis.unmask(0).rename('mask1_defor')
  .addBands(mask2_forest_vs_other_vis.unmask(0).rename('mask2_other'));


Export.image.toDrive({
  image: combinedDisagreementMasks.uint8(),
  description: 'Xingu_Incoherence_Masks_30m',
  folder: 'Article2_Xingu_Panels',
  region: xinguROI,
  scale: 30,
  crs: exportCRS, // FIXED: Replaced mbProj object with explicit clean string
  maxPixels: 1e9
});
*/


//######################################################################
// SECTION 8 : EXPORT INCOHERENCE FRACTION PER 0.02-DEGREE TILE
//######################################################################


// Define our target coarse projection grid (EPSG:4326 at 0.02 degrees)
var coarseProj = ee.Projection('EPSG:4326').atScale(2200);


// Convert our boolean masks to float layers (1 where True, 0 where False)
var mask1_float = mask1_forest_vs_defor_vis.unmask(0).float();
var mask2_float = mask2_forest_vs_other_vis.unmask(0).float();


// Compute the percentage of disagreement per coarse cell
var densityMask1 = mask1_float
  .reduceResolution({ reducer: ee.Reducer.mean(), maxPixels: 65535 })
  .reproject({ crs: coarseProj })
  .multiply(100);


var densityMask2 = mask2_float
  .reduceResolution({ reducer: ee.Reducer.mean(), maxPixels: 65535 })
  .reproject({ crs: coarseProj })
  .multiply(100);


// THE FIX: Create a cookie-cutter mask based on actual MapBiomas valid data
var validDataMask = s2_mb_level1.mask()
  .reduceResolution({
    reducer: ee.Reducer.max(), 
    maxPixels: 65535
  })
  .reproject({ crs: coarseProj });


// Combine bands and apply the true footprint cookie-cutter
var exportDensityImage = densityMask1.rename('pct_incoherence_case1')
  .addBands(densityMask2.rename('pct_incoherence_case2'))
  .updateMask(validDataMask); 


// Get transform tracking from our target projection matrix
var mbTransform = s2_mb_level1.projection().getInfo().transform;
var coarseTransform = [0.01, 0, mbTransform[2], 0, -0.01, mbTransform[5]];


Export.image.toDrive({
  image: exportDensityImage,
  description: 'Amazon_Disagreement_Density_0-01deg',
  folder: 'Article2_Exports',
  region: region, 
  crs: 'EPSG:4326',
  crsTransform: coarseTransform,
  maxPixels: 1e13,
  fileFormat: 'GeoTIFF'
});