Google Earth Engine (GEE) Code to Run and Map the Perimeter-Area Soil Carbon Index (PASCI)

The Perimeter-Area Soil Carbon Index or PASCI could easily be applied to hyperspectral and multispectral images, where wavelengths in the SWIR region are available. PASCI is used to quantify, predict, and map soil organic carbon (SOC). When the PASCI equation is used on a multispectral image such as Sentinel-2 or Landsat, the closest waveband that is available should be applied, while missing wavebands could be skipped.

// This GEE code computes and maps the Perimeter-Area Soil Carbon Index (PASCI), which was developed by Salas and Kumaran, 2023 (see citation below). 

// Salas, E. A. L., & Kumaran, S. S. (2023). Perimeter-Area Soil Carbon Index (PASCI): modeling and estimating soil // organic carbon using relevant explicatory waveband variables in machine learning environment. 
// Geo-Spatial Information Science, 27(6), 1739–1746. https://doi.org/10.1080/10095020.2023.2211612

// Before you run the code:
// 1. Check those lines that are bold. Make sure that those are what you need. Otherwise, change them.
// 2. Your roi should be defined in the imports section of your Earth Engine script
// If roi is not defined elsewhere, uncomment and modify the next line below, and provide your ROI or study area extent.

// var roi = ee.FeatureCollection('your_asset_path_here');  

// Load Sentinel-2 surface reflectance data
var sentinel2 = ee.ImageCollection('COPERNICUS/S2_SR') 
                  .filterBounds(roi)
                  .filterDate('2023-10-01', '2023-12-31') // define your range of date
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20));  // you can change the cloud cover percentage

// Debug: Print the number of images in the collection
print('Number of images in collection:', sentinel2.size());

// Sentinel-2 band information needed
// Landsat images is also possible, but you need to change or edit the specific wavelengths below.
var sentinel2Bands = [
  {name: 'B2', wavelength: 490},  // Blue
  {name: 'B3', wavelength: 560},  // Green
  {name: 'B4', wavelength: 665},  // Red
  {name: 'B5', wavelength: 705},  // Red Edge 1
  {name: 'B6', wavelength: 740},  // Red Edge 2
  {name: 'B7', wavelength: 783},  // Red Edge 3
  {name: 'B8', wavelength: 842},  // NIR
  {name: 'B8A', wavelength: 865}, // NIR narrow
  {name: 'B11', wavelength: 1610},// SWIR1
  {name: 'B12', wavelength: 2190} // SWIR2
];

// Function to calculate PASCI for an image using Earth Engine's band math
var calculatePASCIForImage = function(image) {
  // Select the required bands
  var selectedImage = image.select(sentinel2Bands.map(function(band) { 
    return band.name; 
  }));
  
  // Calculate area using trapezoidal rule
  var area = ee.Image(0);
  
  // Loop to calculate area using the small trapezoids
  for (var i = 0; i < sentinel2Bands.length - 1; i++) {
    var bandWidth = sentinel2Bands[i + 1].wavelength - sentinel2Bands[i].wavelength;
    var bandSum = image.select(sentinel2Bands[i].name).add(image.select(sentinel2Bands[i + 1].name));
    var trapezoidArea = ee.Image(bandWidth).multiply(bandSum).divide(2);
    area = area.add(trapezoidArea);
  }
  
  // Calculate perimeter 
  var perimeter = ee.Image(0);
  
  // Add first reflectance value or vertical line to the perimeter
  var firstReflectance = image.select(sentinel2Bands[0].name);
  perimeter = perimeter.add(firstReflectance);
  
  // Segments along the curve
  for (var i = 0; i < sentinel2Bands.length - 1; i++) {
    var dx = sentinel2Bands[i + 1].wavelength - sentinel2Bands[i].wavelength;
    var dy = image.select(sentinel2Bands[i + 1].name).subtract(image.select(sentinel2Bands[i].name));
    var segmentLength = dy.multiply(dy).add(ee.Image(dx * dx)).sqrt();
    perimeter = perimeter.add(segmentLength);
  }
  
  // Add last reflectance value or vertical lone to the perimeter
  var lastReflectance = image.select(sentinel2Bands[sentinel2Bands.length - 1].name);
  perimeter = perimeter.add(lastReflectance);
  
  // Add the bottom line (total wavelength range) to the perimeter
  var wavelengthRange = sentinel2Bands[sentinel2Bands.length - 1].wavelength - sentinel2Bands[0].wavelength;
  perimeter = perimeter.add(ee.Image(wavelengthRange));
  
  // Calculate P/A ratio for the actual spectrum
  var pOverA = perimeter.divide(area);
  
  // The purpose of this section of the code is to create a baseline or reference spectrum for comparison 
  // with the actual spectrum derived from the Sentinel-2 data. 
  
  // Calculate P/A for simple spectrum (constant reflectance)
  var simpleReflectance = image.select(sentinel2Bands.map(function(band) { return band.name; }))
    .reduce(ee.Reducer.mean());
  
  // Build an image with all bands having the same value (the mean)
  var simpleSpectrum = ee.Image(0);
  for (var i = 0; i < sentinel2Bands.length; i++) {
    simpleSpectrum = simpleSpectrum.addBands(simpleReflectance.rename(sentinel2Bands[i].name));
  }
  simpleSpectrum = simpleSpectrum.select(sentinel2Bands.map(function(band) { return band.name; }));
  
  // Calculate area and perimeter for simple spectrum
  var simpleArea = ee.Image(0);
  for (var i = 0; i < sentinel2Bands.length - 1; i++) {
    var bandWidth = sentinel2Bands[i + 1].wavelength - sentinel2Bands[i].wavelength;
    var bandSum = simpleSpectrum.select(sentinel2Bands[i].name).add(simpleSpectrum.select(sentinel2Bands[i + 1].name));
    var trapezoidArea = ee.Image(bandWidth).multiply(bandSum).divide(2);
    simpleArea = simpleArea.add(trapezoidArea);
  }
  
  var simplePerimeter = simpleReflectance.multiply(2) // First and last vertical lines
    .add(ee.Image(wavelengthRange)); // Add bottom horizontal line
  // Segments along flat line have dx but no dy, so just add wavelength range
  simplePerimeter = simplePerimeter.add(ee.Image(wavelengthRange));
  
  var simplePOverA = simplePerimeter.divide(simpleArea);
  
  // Create a complex zigzag spectrum as maximum reference
  var minReflectance = image.select(sentinel2Bands.map(function(band) { return band.name; }))
    .reduce(ee.Reducer.min());
  var maxReflectance = image.select(sentinel2Bands.map(function(band) { return band.name; }))
    .reduce(ee.Reducer.max());
    
  // Calculate a theoretical maximum P/A for normalization
  var complexPOverA = simplePOverA.multiply(2); // Simplified approach
  
  // Calculate PASCI using the formula from the image
  var pasci = pOverA.subtract(simplePOverA)
    .divide(complexPOverA.subtract(simplePOverA))
    .rename('pasci');
  
  // Clamp values to 0-1 range
  pasci = pasci.clamp(0, 1);
  
  // Add PASCI band to the original image
  return image.addBands(pasci);
};

// Main workflow function
function main() {
  // Apply PASCI calculation to the image collection
  var pasciCollection = sentinel2.map(calculatePASCIForImage);
  
  // Compute a median composite
  var pasciComposite = pasciCollection.select('pasci').median();
  
  // Clip to the roi
  var clippedComposite = pasciComposite.clip(roi);
  
  // Debug: Print the PASCI composite
  print('PASCI Composite:', clippedComposite);
  
  // Visualization parameters
  var visParams = {
    min: 0,
    max: 1,
    palette: ['blue', 'cyan', 'green', 'yellow', 'orange', 'red']
  };
  
  // Add the PASCI composite to the map
  Map.centerObject(roi, 9);
  Map.addLayer(clippedComposite, visParams, 'PASCI Composite');
  
  // Export the result to Drive
  Export.image.toDrive({
    image: clippedComposite,
    description: 'PASCI_Sentinel2_Composite',  // You can change the filename for download
    scale: 10,
    region: roi,
    maxPixels: 1e9
  });
}

// Run the main function
main();

// Create a legend panel
var legend = ui.Panel({
  style: {
    position: 'bottom-right',
    padding: '8px 15px',
    backgroundColor: 'white'
  }
});

// Add a title to the legend
var legendTitle = ui.Label({
  value: 'Soil Carbon (PASCI)',
  style: {
    fontWeight: 'bold',
    fontSize: '14px',
    margin: '0 0 4px 0',
    padding: '0'
  }
});
legend.add(legendTitle);

// Define the color palette and labels. If you want, change the labels and range of values.
var palette = ['blue', 'cyan', 'green', 'yellow', 'orange', 'red'];
var labels = [
  'Very Low SOC (0.0 - 0.2)',
  'Low SOC (0.2 - 0.4)',
  'Moderate SOC (0.4 - 0.6)',
  'High SOC (0.6 - 0.8)',
  'Very High SOC (0.8 - 1.0)',
  'Extremely High SOC (1.0)'
];

// Add color boxes and labels to the legend
for (var i = 0; i < palette.length; i++) {
  var colorBox = ui.Label({
    style: {
      backgroundColor: palette[i],
      padding: '8px',
      margin: '0 0 4px 0'
    }
  });

  var label = ui.Label({
    value: labels[i],
    style: {
      margin: '0 0 4px 8px'
    }
  });

  legend.add(colorBox).add(label);
}

// Add the legend to the map
Map.add(legend);