Harmonized Landsat 5–9 NDVI Trend (1984–today, Mann-Kendall + Sen Slope)

// ============================================================
// Harmonized Landsat NDVI Trend (1984 -> today)
// Mann-Kendall + Sen slope per pixel, p<0.05 significance
// ============================================================
var startDate = '1984-01-01';
var endDate   = ee.Date(Date.now());

// --- Cloud / shadow mask via QA_PIXEL bitmask
function maskL2(img) {
  var qa = img.select('QA_PIXEL');
  var dilatedCloud = 1 << 1, cloud = 1 << 3, shadow = 1 << 4;
  var mask = qa.bitwiseAnd(dilatedCloud).eq(0)
              .and(qa.bitwiseAnd(cloud).eq(0))
              .and(qa.bitwiseAnd(shadow).eq(0));
  return img.updateMask(mask)
            .multiply(0.0000275).add(-0.2)
            .copyProperties(img, ['system:time_start']);
}

// --- Harmonize each sensor to a common (red, nir) pair
function harmonize(coll, redB, nirB) {
  return coll.map(maskL2).map(function(i){
    return i.normalizedDifference([nirB, redB]).rename('NDVI')
            .copyProperties(i, ['system:time_start']);
  });
}
var L5 = harmonize(ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').filterBounds(aoi), 'SR_B3', 'SR_B4');
var L7 = harmonize(ee.ImageCollection('LANDSAT/LE07/C02/T1_L2').filterBounds(aoi), 'SR_B3', 'SR_B4');
var L8 = harmonize(ee.ImageCollection('LANDSAT/LC08/C02/T1_L2').filterBounds(aoi), 'SR_B4', 'SR_B5');
var L9 = harmonize(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2').filterBounds(aoi), 'SR_B4', 'SR_B5');
var ndviAll = L5.merge(L7).merge(L8).merge(L9).filterDate(startDate, endDate);

// --- Annual maxima (anti-cloud composite)
var years = ee.List.sequence(1984, ee.Date(endDate).get('year'));
var annual = ee.ImageCollection(years.map(function(y){
  y = ee.Number(y);
  var img = ndviAll
    .filter(ee.Filter.calendarRange(y, y, 'year'))
    .max().rename('NDVI');
  return img.set('year', y).set('system:time_start', ee.Date.fromYMD(y,7,1).millis());
}));

// --- Mann-Kendall tau + Sen slope per pixel
var afterFilter = ee.Filter.lessThan({leftField:'system:time_start', rightField:'system:time_start'});
var joined = ee.ImageCollection(ee.Join.saveAll('after').apply({
  primary: annual, secondary: annual, condition: afterFilter
}));
var sign = ee.ImageCollection(joined.map(function(curr){
  var afters = ee.ImageCollection.fromImages(curr.get('after'));
  return afters.map(function(a){
    return ee.Image(a).neq(curr).multiply(ee.Image(a).subtract(curr).where(a.eq(curr), 0).divide(ee.Image(a).subtract(curr).abs()));
  });
}).flatten());
var kendall = sign.reduce('sum', 2).rename('kendall');
var slope = annual.reduce(ee.Reducer.sensSlope()).select('slope').rename('sen_slope');

// --- Significance (n=count, var = n(n-1)(2n+5)/18)
var n = annual.count();
var varS = n.multiply(n.subtract(1)).multiply(n.multiply(2).add(5)).divide(18);
var z = kendall.divide(varS.sqrt());
var pSig = z.abs().gt(1.96);   // ~p<0.05 two-tailed

// --- Visualisation
Map.centerObject(aoi, 7);
Map.addLayer(slope.updateMask(pSig), {min:-0.01, max:0.01, palette:['#dc2626','#f5f5f5','#16a34a']}, 'NDVI Sen slope (sig)');
Map.addLayer(pSig.updateMask(pSig), {palette:['#000000']}, 'p<0.05 mask', false);

// --- Charts + zonal stats
print(ui.Chart.image.series({
  imageCollection: annual, region: aoi, reducer: ee.Reducer.mean(),
  scale: 1000, xProperty: 'system:time_start'
}).setOptions({title: 'Annual NDVI maxima (AOI mean)', vAxis: {title: 'NDVI'}}));

var zonal = slope.addBands(pSig.rename('sig')).reduceRegions({
  collection: aoi.geometry().bounds().coveringGrid('EPSG:4326', 10000),
  reducer: ee.Reducer.mean(), scale: 1000, tileScale: 4
});

// --- Exports
Export.image.toDrive({image: slope.clip(aoi), description: 'NDVI_SenSlope',
  folder: 'gee_exports', region: aoi, scale: 30, maxPixels: 1e13,
  fileFormat: 'GeoTIFF', formatOptions: {cloudOptimized: true}});
Export.image.toDrive({image: pSig.clip(aoi),  description: 'NDVI_pSig',
  folder: 'gee_exports', region: aoi, scale: 30, maxPixels: 1e13,
  fileFormat: 'GeoTIFF', formatOptions: {cloudOptimized: true}});
Export.table.toDrive({collection: zonal, description: 'NDVI_Trend_Zonal_CSV',
  folder: 'gee_exports', fileFormat: 'CSV'});

# Python equivalent (geemap). See JS for full inline docs.
import ee, geemap
ee.Initialize()
m = geemap.Map()
# ... build harmonized collection as in the JS version, then:
# slope = annual.reduce(ee.Reducer.sensSlope()).select('slope')
# m.addLayer(slope, {'min': -0.01, 'max': 0.01,
#   'palette': ['#dc2626','#f5f5f5','#16a34a']}, 'NDVI Sen slope')