【GEE】遥感数据趋势分析Sen+mk

Map.centerObject(table);

// 定义时间范围
var stary = 2001, endy = 2023;
//NDVI图像集合
var NDVICL = ee.ImageCollection(ee.List.sequence(stary, endy).map(function(year) {
  // 定义每年的开始和结束日期
  var startd = ee.Date.fromYMD(year, 1, 1);
  var endd = ee.Date.fromYMD(year, 12, 31);
  // 加载NDVI数据集，按日期过滤，并计算最大值
  return ee.ImageCollection('MODIS/006/MOD13A1')
    .filterDate(startd, endd)
    .select('NDVI')
    .max()
    .addBands(ee.Image.constant(year).toFloat().rename('year'));
}));
print(NDVICL)
// 计算Sen斜率并剪裁
var senSlope = NDVICL.select(['NDVI', 'year'])
    .reduce(ee.Reducer.sensSlope())
    .clip(table);
// 可视化Sen斜率
var vis = {bands: ['slope'], min: -0.5, max: 0.5, palette: ['red', 'white', 'green']};
Map.addLayer(senSlope.select('slope'), vis, "Sen斜率 ");

// MK检验 
// create ones image
var ones = ee.Image(1);
// create zeros image
var zeros = ee.Image(0);
var eps = 0.01; // Set a precision value
// 创建一个列表，包含年度NDVI图像的list
var listofimg = NDVICL.toList(NDVICL.size());
// 计算MK检验趋势期望均值
var mkm = ee.ImageCollection(ee.List.sequence(0, listofimg.size().subtract(2)).map(
    function (xi) {
        return ee.ImageCollection.fromImages(
            ee.List.sequence(ee.Number(xi).add(1), listofimg.size().subtract(1)).map(
                function (xj) {
                    var diff_ij = ee.Image(listofimg.get(xj)).select('NDVI')
                        .subtract(ee.Image(listofimg.get(xi)).select('NDVI'));
                    var diff = diff_ij.where(diff_ij.abs().lt(eps), zeros);
                    return diff.where(diff.gt(0), ones)
                        .where(diff.lt(0), ones.multiply(-1)).rename('ES');
                     })).sum();
    })).sum();
// 时序长度 
var n = ee.ImageCollection(NDVICL)
    .reduce('count')
    .select('NDVI_count')
    .rename('count');
// 计算方差 var(S) = (n * (n - 1) * (2 * n + 5)) / 18
var varS = n.multiply(
    n.subtract(1).multiply(
        n.multiply(2).add(5)))
    .divide(18).rename('varS');

// 估计z统计量
var mkz = mkm.where(mkm.abs().lt(eps), 0);
var zcore = mkz
    .where(mkz.gt(eps), mkz.subtract(1).divide(varS.sqrt()))
    .where(mkz.lt(-eps), mkz.add(1).divide(varS.sqrt())).rename('zcore')
    .clip(table);

// 100(1 - alpha)% 双边置信区间的Sen's Slope
//alpha = 0.05，则Z[1 - alpha / 2] = Z[0.975] = 1.960
// 计算Sen斜率的显著性
var signif = senSlope.select(0).lt(0).and(zcore.abs().gte(1.960)).rename('signif');
var siginc = senSlope.select(0).gt(0).and(zcore.abs().gte(1.960)).rename('increase');

signif.updateMask(signif);
siginc.updateMask(siginc);

signif = signif.multiply(-1).add(siginc).clip(table);

// 可视化MK Z统计量
Map.addLayer(zcore, {min: -3, max: 3, palette: ['blue', 'white', 'red']}, 'Z值 ');

// 可视化显著的趋势变化
Map.addLayer(signif, {min: -1, max: 1, palette: ['#FFA07A', 'white', '#3CB371']}, '显著趋势 ');