吴敦义党主席之路面临严峻挑战

Question

I'm trying to run a Mann-Kendall analysis on spatial NDVI timeseries for a number of sites. As part of the analyses, I have been using DynamicWorld to mask out undesirable areas (e.g. snow and ice, bare ground, water, etc.). I was initially following a GEE forum tutorial for the MK analysis, but myself and others have found issues with the way the sum of the signs were calculated. I tried to come up with a solution myself:

// Function to calculate sign between two images
// If NDVI at time j is more than NDVI at time i, the sign is 1
    // if opposite is true, sign is -1
    // if equal, sign is 0
var sign = function(i, j) {
  var diff = ee.Image(j).subtract(i);
  // Applying DynamicWorld masking to starting image
  var zero = ee.Image(0).updateMask(classmask.not());
  return zero
    .where(diff.gt(0), 1)
    .where(diff.lt(0), -1)
    .int();
};

But this has created issues where pixels previously masked by DynamicWorld are exported as 0 values within the CSV, despite appearing as masked when added to the map before export. I've tried a number of solutions, including stating the mask again before export or setting as nodata values in export, but with no luck. This seems to be an issue specifically with Mann-Kendall S statistic calculation, because there is no pixel inflation in later stages when exporting p values or Sen's Slope.

Full code below.

// Load site geometry
var allsites = ee.FeatureCollection('users/USERNAME/allsites');

// Select one site (example)
var siteId = 'Border_Meuse';
var site = allsites.filter(ee.Filter.eq('site_id', siteId));
var geometry = site.geometry().dissolve();

///////////////////////
//                   //
//    DynamicWorld   //
//                   //
///////////////////////

var dworld = ee.ImageCollection("GOOGLE/DYNAMICWORLD/V1")
.filterDate('2025-08-06', '2025-08-06')
.filterBounds(site);

var allclass = ['water','trees','grass','flooded_vegetation','crops','shrub_and_scrub',
  'built','bare','snow_and_ice'];
var undesirableclass = ['water', 'crops', 'built', 'bare', 'snow_and_ice'];

var siteCompositor = function(imageCollection) {
  return imageCollection
    .select(allclass)
    .median()
    .copyProperties(imageCollection.first(), imageCollection.first().propertyNames());
};

var dworldcomp = ee.Image(siteCompositor(dworld)).clip(site.geometry());

var classprob = dworldcomp.reduce(ee.Reducer.sum());

var dworldsites = dworldcomp.divide(classprob);

// Sum the probabilities of the undesirable classes
var classprobsum = dworldsites.select(undesirableclass).reduce(ee.Reducer.sum());

// Mask where the total probability > 0.75
var classmask = classprobsum.gt(0.75);

///////////////////////////////////
//                               //
//        BUILDING IMAGES        //
//                               //
///////////////////////////////////

// Define the year range
var startYear = ee.Number(site.first().get('start_year'));
var endYear = 2024;
var years = ee.List.sequence(startYear, endYear);

var yearList = years.getInfo().filter(function(y) {
  return y !== 2003;
});

// Build image list
var images = yearList.map(function(y) {
  var path = 'users/AidanHoutenUCL/ndvi_rasters/' + siteId + '/' + y;
  var image = ee.Image(path).select('NDVI')
    .set('year', y)
    .set('system:time_start', ee.Date.fromYMD(y, 6, 1).millis());
  return image;
});

// Convert to ee.ImageCollection
var coll = ee.ImageCollection(images);

var cleanedNDVI = coll.map(function(image) {
  return image.updateMask(classmask.not());
});

////////////////////////////////
//                            //
//        MANN-KENDALL        //
//                            //
////////////////////////////////

// MK analysis is used to detect monotonic trends i.e. a consistent increase/decrease over time

// Define temporal filter to pass all the chronologically later images
var afterFilter = ee.Filter.lessThan({
  leftField: 'year',
  rightField: 'year'
});

// Join the collection to itself based on year
// each image will store all the images that come after it in an 'after' property
var joined = ee.ImageCollection(ee.Join.saveAll('after').apply({
  primary: cleanedNDVI,
  secondary: cleanedNDVI,
  condition: afterFilter
}));

// Function to calculate sign between two images
// If NDVI at time j is more than NDVI at time i, the sign is 1
    // if opposite is true, sign is -1
    // if equal, sign is 0
var sign = function(i, j) {
  var diff = ee.Image(j).subtract(i);
  // Applying DynamicWorld masking to starting image
  var zero = ee.Image(0).updateMask(classmask.not());
  return zero
    .where(diff.gt(0), 1)
    .where(diff.lt(0), -1)
    .int();
};

// Calculate Mann-Kendall S Statistic
var kendall = ee.ImageCollection(joined.map(function(current) {
  var afterCollection = ee.ImageCollection.fromImages(current.get('after'));
  return afterCollection.map(function(image) {
    return ee.Image(sign(current, image));
  });
}).flatten()).reduce('sum', 2);

var kendall=kendall.clip(geometry)

// Visualize Mann-Kendall S Statistic
Map.centerObject(geometry);
Map.addLayer(kendall, {min: -100, max: 100, palette: ['red', 'white', 'green']}, 'Mann-Kendall S');

// Export result
Export.image.toDrive({
  image: kendall,
  description: 'MannKendall_' + siteId,
  folder: 'GEE_MK',
  fileNamePrefix: 'MannKendall_' + siteId,
  region: geometry,
  scale: 30,
  maxPixels: 1e13
});

////////////////////////////////
//                            //
//    VARIANCE CALCULATION    //
//                            //
////////////////////////////////

// Detects tied NDVI values
var groups = cleanedNDVI.map(function(i) {
  var matches = cleanedNDVI.map(function(j) {
    return i.eq(j);
  }).sum();
  return i.multiply(matches.gt(1));
});

// Compute tie group sizes in a sequence.  The first group is discarded.
var group = function(array) {
  var length = array.arrayLength(0);
  // Array of indices.  These are 1-indexed.
  var indices = ee.Image([1])
      .arrayRepeat(0, length)
      .arrayAccum(0, ee.Reducer.sum())
      .toArray(1);
  var sorted = array.arraySort();
  var left = sorted.arraySlice(0, 1);
  var right = sorted.arraySlice(0, 0, -1);
  // Indices of the end of runs.
  var mask = left.neq(right)
  // Always keep the last index, the end of the sequence.
      .arrayCat(ee.Image(ee.Array([[1]])), 0);
  var runIndices = indices.arrayMask(mask);
  // Subtract the indices to get run lengths.
  var groupSizes = runIndices.arraySlice(0, 1)
      .subtract(runIndices.arraySlice(0, 0, -1));
  return groupSizes;
};

// See equation 2.6 in Sen (1968) - implements a correction factor for ties in the data
// The sum of the below over all tied groups then subtracted from the general variance formula
var factors = function(image) {
  return image.expression('b() * (b() - 1) * (b() * 2 + 5)');
};

// Identifies runs of tied values (e.g. NDVI = 0.4 for 3 years in a row)
// Computes the length of each run
var groupSizes = group(groups.toArray());
var groupFactors = factors(groupSizes);
var groupFactorSum = groupFactors.arrayReduce('sum', [0])
      .arrayGet([0, 0]);

// Per-pixel Kendall's variance
var count = joined.count();
var kendallVariance = factors(count)
    .subtract(groupFactorSum)
    .divide(18)
    .float();

// Map.addLayer(kendallVariance, {min: 0, max: 10000}, 'kendallVariance');

////////////////////////////////////
//                                //
//    SIGNIFICANCE CALCULATION    //
//                                //
////////////////////////////////////

// Compute Z-statistics per-pixel based on the sign of the trend and its variance
var zero = kendall.multiply(kendall.eq(0));
var pos = kendall.multiply(kendall.gt(0)).subtract(1);
var neg = kendall.multiply(kendall.lt(0)).add(1);
var z = zero
    .add(pos.divide(kendallVariance.sqrt()))
    .add(neg.divide(kendallVariance.sqrt()));

// Map.addLayer(z, {min: -2, max: 2}, 'z');

// Function to convert Z to p-value using a cumulative distribution function
function eeCdf(z) {
  return ee.Image(0.5)
      .multiply(ee.Image(1).add(ee.Image(z).divide(ee.Image(2).sqrt()).erf()));
}

// Inverse cumulative distribution function
// Calculates the Z-score corresponding to a given p-value
function invCdf(p) {
  return ee.Image(2).sqrt()
      .multiply(ee.Image(p).multiply(2).subtract(1).erfInv());
}

// Compute P-values
var p = ee.Image(1).subtract(eeCdf(z.abs()));
Map.addLayer(p, {min: 0, max: 1}, 'p');

// Pixels that can have the null hypothesis (there is no trend) rejected
Map.addLayer(p.lte(0.025), {min: 0, max: 1}, 'significant trends');

// Export
Export.image.toDrive({
  image: p,
  description: 'PValues_' + siteId,
  folder: 'GEE_P',
  fileNamePrefix: 'PValues_' + siteId,
  region: geometry,
  scale: 30,
  maxPixels: 1e13
});

///////////////////////////////
//                           //
//        SEN'S SLOPE        //
//                           //
///////////////////////////////

// Define slope function between two images (NDVI diff / time diff in years)
var slope = function(i, j) {
  return ee.Image(j).subtract(i)
      .divide(ee.Image(j).date().difference(ee.Image(i).date(), 'year'))  // Use years
      .rename('slope')
      .float();
};

// Apply slope function to all image pairs
var slopeImages = ee.ImageCollection(joined.map(function(current) {
  var afterCollection = ee.ImageCollection.fromImages(current.get('after'));
  return afterCollection.map(function(image) {
      return ee.Image(slope(current, image));
  });
}).flatten());

// Reduce to median slope (Sen's slope)
var sensSlope = slopeImages.reduce(ee.Reducer.median()).rename('Sen_Slope');
var sensSlopeClipped = sensSlope.clip(geometry);

// Add to map
Map.addLayer(sensSlopeClipped, {min: -0.01, max: 0.01, palette: ['red', 'white', 'green']}, 'Sen\'s Slope');

// Export
Export.image.toDrive({
  image: sensSlopeClipped,
  description: 'SenSlope_' + siteId,
  folder: 'GEE_SS',
  fileNamePrefix: 'SenSlope_' + siteId,
  region: geometry,
  scale: 30,
  maxPixels: 1e13
});