基于GEE的bfastmonitor的改编

bfast算法简介:

bfast是季节和趋势分量的迭代中断检测系列。季节性中断是一个函数，它将时间序列的迭代分解为趋势分量、季节分量和剩余分量，并在分解后的时间序列分量中检测到显著中断。

简单来说,就是分析一组数据,检测他们出现突变的时间点及其程度。

bfast的r语言包：https://r-forge.r-project.org/projects/bfast/

在gee上的bfastmonitor能很好地对遥感影像进行计算，其APP：BFASTmonitorGEE

该算法在github上的monitor代码，其包含了遥感影像的处理、NDWI等指数的计算、再利用bfast算法对其分析：https://github.com/bfast2/geeBfastMonitor

但是博主需要对gee上已有数据集进行分析，于是博主删去bfastMonitor上的指数计算部分，再掩膜掉数据集中的缺失值，使得该算法能直接对数据集进行计算，

博主基于bfastMonitor改编的算法：https://code.earthengine.google.com/a31675794094ab05c6bbaa52b7f12cbcz

使用步骤：

1.选择自己感兴趣的区域。

2.选择数据，这个数据一定得是imageCollection数据集，并且有足够的影像。

3.选择想要计算的波段，该算法一次性只能计算一种波段。

4.保存并导出图像。

全代码，总共有五百多行：

//1.这里上次自己感兴趣的区域和数据集
var roi=ee.FeatureCollection('users/20220404g/xizang')var imageCollection2=ee.ImageCollection("LANDSAT/LC8_L1T_32DAY_NDWI")//done
var imageCollection7=ee.ImageCollection("ECMWF/ERA5_LAND/MONTHLY")//debug done
var imageCollection8=ee.ImageCollection("NASA/GRACE/MASS_GRIDS/MASCON_CRI")//doneprint(imageCollection2)
print(imageCollection7)//parameter you could change，2.这里需要修改为自己选的数据集，包括下面的参数，可以根据自己的需要修改
var collection=imageCollection7.filterBounds(roi)//you must change
var historyStart = '2011-01-01'
var historyEnd ='2014-12-31'
var monitoringStart ="2015-01-01"
var monitoringEnd ="2016-12-31"
var yearStart=2011
var yearEnd=2017
var dependent = 'snow_cover';var h =0.25
var hxh = h
var period = 10
var alpha =0.05
var atAlpha = 1 - alpha;
var magnitudeThreshold =-0.000000000000000000000000000000121
var harmonics =1//better not to change，这下面的尽量不要碰
var roio = roi
var selectBand=function(image){return image.select(dependent)
}var collectionSeBand=collection.map(selectBand);var histCollection = collectionSeBand.filterDate(historyStart, historyEnd);
var moniCollection  = collectionSeBand.filterDate(monitoringStart, monitoringEnd);//Check the h parameter
if (h == 0.25) {h = 0;
} else if (h == 0.5){h = 1;
}else if (h == 1){h = 2;
}else {alert("Error: h parameter can either be 0.25, 0.5, or 1");
}///Check alpha. Alpha be equal or less than 0.05if (alpha > 0.05) {alert("Error: alpha parameter set too large, try alpha = 0.05 or less");
}var tTable = [0.950, 0.951, 0.952, 0.953, 0.954, 0.955, 0.956, 0.957, 0.958, 0.959, 0.960, 0.961, 0.962, 0.963, 0.964, 0.965, 0.966, 0.967, 0.968, 0.969,
0.970, 0.971 ,0.972, 0.973, 0.974, 0.975, 0.976, 0.977, 0.978, 0.979, 0.980, 0.981, 0.982, 0.983, 0.984, 0.985, 0.986, 0.987, 0.988, 0.989, 0.990 ,0.991 ,0.992,
0.993 ,0.994, 0.995, 0.996, 0.997, 0.998, 0.999];var alphaIndex = tTable.indexOf(atAlpha)
if (alphaIndex < 0 || alphaIndex > 49) {alert("Error: for critical values, we only have for alpha parameter ranging from 0.001 to 0.05, try alpha from that range");
}if (period == 2) {period  = 0;
} else if (period  == 4){period  = 1;
}else if (period  == 6){period  = 2;
}else if (period  == 8){period  = 3;
}else if (period  == 10){period  = 4;
}else {alert("Error: for period parameter, we only have 2, 4, 6, 8,10. Choose one of these values");
}
///
// DEFINE KEY FUNCTIONS FOR PRE=PROCESSING
///// Make a list of harmonic frequencies to model.
// These also serve as band name suffixes.
var harmonicFrequencies = ee.List.sequence(1, harmonics);// Function to get a sequence of band names for harmonic terms.
var constructBandNames = function(base, list) {return ee.List(list).map(function(i) {return ee.String(base).cat(ee.Number(i).int());});
};// Construct lists of names for the harmonic terms.
var cosNames = constructBandNames('cos_', harmonicFrequencies);
var sinNames = constructBandNames('sin_', harmonicFrequencies);// Independent variables.
var independents = ee.List(['constant']).cat(cosNames).cat(sinNames);// Functions to add a time band.//I add a transform
var addDependents = function(image) {// Compute time in fractional years since the epoch.var years = image.date().difference('1970-01-01', 'year');var timeRadians = ee.Image(years.multiply(2 * Math.PI)).rename('t');var constant = ee.Image(1);return image.addBands(constant).addBands(timeRadians.float());
};// Function that returns the year (with a fractional part to indicate progression within
// a year) that corresponds to the given unix time.
function unixToYear(unixTime) {return ee.Number(unixTime).divide(365.25 * 24 * 3600 * 1000).add(1970);
}// Function to compute the specified number of harmonics
// and add them as bands.  Assumes the time band is present.
var addHarmonics = function(freqs) {return function(image) {// Make an image of frequencies.var frequencies = ee.Image.constant(freqs);// This band should represent time in radians.var time = ee.Image(image).select('t');// Get the cosine terms.var cosines = time.multiply(frequencies).cos().rename(cosNames);// Get the sin terms.var sines = time.multiply(frequencies).sin().rename(sinNames);return image.addBands(cosines).addBands(sines);};
};
//print("done")/HISTORY PERIOD var harmonicLandsat = histCollection.map(addDependents).map(addHarmonics(harmonicFrequencies))// The output of the regression reduction is a 4x1 array image.
var harmonicTrend = harmonicLandsat.select(independents.add(dependent)).reduce(ee.Reducer.linearRegression(independents.length(), 1));// Turn the array image into a multi-band image of coefficients.
var harmonicTrendCoefficients = harmonicTrend.select('coefficients').arrayProject([0]).arrayFlatten([independents]);// 3. CALCULATE COMMULATIVE SUM OF RESIDUALS// Compute fitted values.//add 5 band
var fittedHarmonic = harmonicLandsat.map(function(image) {//var transform=image.projection()return image.addBands(image.select(independents).multiply(harmonicTrendCoefficients).reduce('sum').rename('fitted'))//.reproject({crs:image.projection()});
});
//print(fittedHarmonic)
// Compute the residuals
var computResiduals = function(image){return image.addBands(image.select(dependent).subtract(image.select('fitted')).reduce('sum').rename('residual'));}var residuals = fittedHarmonic.map(computResiduals);// Compute squared residuals
var computsquaredResiduals = function(image){//var transform=image.projection()return image.addBands(image.select('residual').pow(2).rename('squaredResidual'))//.reproject({crs:image.projection()});
}
//add 7 band
var squaredResiduals = residuals.map(computsquaredResiduals);
//print(squaredResiduals)
//Map.addLayer(squaredResiduals.select(dependent), 'value of test band');
var rndmi = residuals.select(dependent)
var nmdi = rndmi.toArray();
//calculate sigma i.e. sigma = sqrt(sum(e^2)/fm$df.residual)
var residu = residuals.select('residual')
var res = residu.toArray();
var resid = squaredResiduals.select('squaredResidual')
var residarray = resid.toArray();
var imageAxis = 0;
var ResidSum = residarray.arrayReduce("sum", [imageAxis]);
var coefArray = harmonicTrendCoefficients.toArray();
var coefLength = coefArray.arrayLength(imageAxis);
var histtimeseiesLength = residarray.arrayLength(imageAxis);
var dfresiduals = histtimeseiesLength.subtract(coefLength);
var ResidSumnorma = ResidSum.divide(dfresiduals);
var sigma = ResidSumnorma.sqrt();/// calculate K, which is N obs in historical sample multiplied by h parameter
var himage = ee.Image(hxh);
var ksize =  himage.multiply(histtimeseiesLength).floor();// Calculate mosum of residuals for history period /// Concatenate two images: zer0 and residuals
var zer0 = ee.Image([0]).toArray().float();
var bandAxis = 1
var fi = zer0.arrayRepeat(bandAxis, 1);
var catImage = ee.ImageCollection([fi,res]);
var concat = catImage.toArrayPerBand(0);
var ImToArr = concat.toArray();/// ...now calculate Cumsum
var histcumsum = ImToArr.arrayAccum(imageAxis, 'sum');
var histimZer0 = ee.Image(0).int();
var histdubImZer0 = histimZer0.arrayRepeat(imageAxis,ksize.int());
var histfx10 = histdubImZer0.arrayRepeat(bandAxis, 1);
var histverc= histfx10.toArray().float();
var histshiftedCumR = ee.ImageCollection([histverc, histcumsum]);
var histShfCumSumRight1 = histshiftedCumR.toArrayPerBand(imageAxis);
var histShfCumSumRight = histShfCumSumRight1.toArray()
var histshiftedCumL = ee.ImageCollection([histcumsum, histverc]);
var histShfCumSumLeft = histshiftedCumL.toArrayPerBand(imageAxis);
var histShiftDifference = histShfCumSumLeft.subtract(histShfCumSumRight);
var histmsk = histShiftDifference.arrayMask(histShfCumSumLeft);
var histter = ee.Image(0).float().toArray();
var histtesser = histter.arrayRepeat(bandAxis, 1);
var histtessermaske = ee.ImageCollection([histtesser, histmsk]);
var histtessermsk = histtessermaske.toArrayPerBand(imageAxis);// ...now slash off unwanted first part of Cumsum of residuals to get MOSUM residuals
var histcotx2 = histtessermsk.arrayLength(imageAxis);
var histcumKsizeDiff = histcotx2.subtract(ksize).int();
var histknsiz1 = histimZer0.arrayRepeat(imageAxis, ksize.int());
var histknsiz2 = histknsiz1.arrayRepeat(bandAxis, 1);
var histknsiz =histknsiz2.toArray();var histdcumn = ee.Image(1).multiply(histcumKsizeDiff);
var histimOnes = ee.Image(1).int();
var histknm1 = histimOnes.arrayRepeat(imageAxis, histdcumn);
var histknm2 = histknm1.arrayRepeat(bandAxis, 1);
var histknm = histknm2.toArray();var histshmasker = ee.ImageCollection([histknsiz, histknm]);
var histShfmask = histshmasker.toArrayPerBand(imageAxis);
var histResmosum = histtessermsk.arrayMask(histShfmask);/// ...standardise
var sigmaStanda = sigma.multiply(histtimeseiesLength.sqrt());var duSigmaStanda = sigmaStanda.arrayRepeat(bandAxis,1);
var bleng = histResmosum.arrayLength(imageAxis);
var xSigmaStanda = duSigmaStanda.arrayRepeat(imageAxis,bleng).toArray();var histResmosumStanda = histResmosum.divide(xSigmaStanda);
//print("done2")///CRITICAL VALUE TABLE ///var criticalVTable = [[[1.227627,1.336231,1.341087,1.341657,1.341825],[1.687323,1.886331,1.899584,1.901299,1.902003],[2.224088,2.704437,2.737148,2.742879,2.745928]],
[[1.230670,1.338791,1.343916,1.344138,1.344391],[1.691601,1.890238,1.903723,1.905166,1.905759],[2.231672,2.713164,2.743205, 2.749723,2.753326]],
[[1.232765,1.341468,1.346242,1.346456,1.346603],[1.696334,1.895051,1.907863,1.909372,1.910032],[2.238556,2.722308,2.750227,2.757492,2.760331]],[[1.235641,1.344179,1.348399,1.348852,1.349151],
[1.701584,1.899687,1.912100,1.913716,1.914301],[2.246716,2.730136,2.757795,2.765451,2.767957]],[[1.238478,1.346586,1.351096,1.351554,1.351786],[1.705517,1.903961,1.916485,1.917941,1.918521],
[2.254955,2.737834,2.766094,2.772398,2.774493]],[[1.241981,1.349207,1.353765,1.353986,1.354179],[1.711073,1.908307, 1.920913,1.922321,1.923639],[ 2.263672,2.745290,2.772925,2.780656,2.783772]],
[[1.244816,1.352020,1.356149,1.356345,1.356684],[1.716266,1.912950,1.926254,1.927599,1.928130],[2.271981,2.753501,2.782581,2.788330,2.790409]],[[1.248790,1.354620,1.358915,1.359235,1.359487],
[1.720193,1.917516,1.931100,1.932476,1.933184],[2.280294,2.761655,2.789957,2.795878,2.797913]],[[1.252395,1.357197,1.361947,1.362265,1.362569],[1.725837,1.922129,1.936171,1.937710,1.938192],
[2.289980,2.769950,2.797966,2.804613,2.808125]],[[1.254989,1.360500,1.365236,1.365600,1.365772],[1.730801,1.927967,1.941870,1.943303,1.943724],[2.301392,2.777844,2.808374,2.813311,2.815859]],
[[1.258229,1.363725,1.368277,1.368505,1.368863],[1.736322,1.933243,1.947243,1.948883,1.949207],[2.310398,2.788013,2.816225,2.821322,2.824270]],[[1.262239,1.366679,1.371643,1.372149,1.372374],
[1.741964,1.939180,1.952111,1.953539,1.954109],[2.320316,2.796900,2.825448,2.831913,2.834508]],[[1.265966,1.370524,1.374604,1.374777,1.374852],[1.747394,1.945610,1.957768,1.959244,1.959426],
[2.329219,2.807824,2.835817,2.840917,2.843434]],[[1.269240,1.373742,1.377815,1.378134,1.378315],[1.753500,1.950991,1.963491,1.964871,1.965069],[2.339306,2.817101,2.846078,2.851046,2.853604]],
[[1.272641,1.376577,1.381153,1.381548,1.381751],[1.758805,1.957234,1.969616,1.970576,1.970974],[2.350632,2.828103,2.856533,2.860616,2.862433]],[[1.276040,1.380221,1.384750,1.385074,1.385378],
[1.765472,1.963210,1.974486,1.975609,1.975930],[2.362393,2.839029,2.866151,2.871058,2.872654]],[[1.279592,1.383943,1.388012,1.388368,1.388473],[1.772518,1.969691,1.980481,1.981478,1.981607],
[2.373682,2.849856,2.875078,2.878643,2.880942]],[[1.284859,1.387693,1.390913,1.391173,1.391456],[1.779402,1.975368,1.985724,1.987007,1.987355],[2.383993,2.859437,2.885461,2.889821,2.891402]],
[[1.289184,1.390484,1.394646,1.395081,1.395330],[1.788029,1.981488,1.992285,1.993116,1.993443],[2.396559,2.872062,2.896479,2.900194,2.901336]],[[1.293119,1.394082,1.398722,1.398860,1.399112],
[1.795241,1.987604,1.997985, 1.998916,1.999079],[2.409003,2.881957,2.906464,2.911006,2.912487]],[[1.297743,1.398436,1.402654,1.403074,1.403188],[1.802223,1.994707,2.004819,2.006503,2.006985],
[ 2.420513,2.894229,2.918304,2.920887,2.922340]],[[1.302930,1.402841,1.406762,1.407276,1.407490],[1.809158,2.001483,2.012325,2.013257,2.013485],[2.431878,2.905888,2.927930,2.931367,2.933102]],
[[1.307095,1.407343,1.411457,1.411639,1.411814],[1.816827,2.009906,2.019668,2.020551,2.020959],[2.442134,2.918303,2.939940,2.942373,2.943662]],[[1.311586,1.411896,1.415374,1.415582,1.415698],
[1.824857,2.018009,2.027714,2.028695, 2.029367],[2.455465,2.928892,2.951336,2.955135,2.955942]],[[1.317376,1.415971,1.419370,1.419562,1.419777],[1.833343,2.026390,2.035115,2.035950,2.036448],
[2.468240,2.941650,2.962908,2.965146,2.966898]],[[1.323352,1.420220,1.423625,1.423804,1.423819],[1.841864,2.034022,2.042662,2.044230,2.044388],[2.483054,2.955380,2.976538,2.979340, 2.980014]],
[[1.327864,1.425241,1.428682,1.428849,1.428957],[1.851771,2.042079,2.052127,2.053173,2.053381],[2.500143,2.968288,2.989569,2.992682,2.994808]],[[1.333507,1.430859, 1.433526,1.433628,1.433639],
[1.861211,2.052014,2.058428,2.059145,2.059377],[2.512775,2.983458,3.002817,3.007144,3.008677]],[[1.339192,1.435554,1.438252,1.438392,1.438405],[1.869505,2.058946,2.065555,2.066094,2.066162],
[2.527290,3.000886,3.017675,3.021515,3.022463]],[[1.344926,1.440531,1.442914,1.442974, 1.443076],[1.878587,2.066095,2.074048,2.074708,2.074738],[2.544121,3.014506,3.031640,3.033255,3.033942]],
[[1.350273,1.445440,1.448123,1.448228,1.448236],[1.886953,2.075533,2.082448,2.082848,2.082870],[2.559929,3.030217,3.045733,3.048442,3.049289]],[[1.356197,1.450892,1.453143,1.453262,1.453311],
[1.899271,2.085649,2.092158,2.092533,2.092569],[2.579905,3.045143,3.062709,3.065217,3.065598]],[[1.362590,1.456028,1.458898,1.458992, 1.459029],[1.909990,2.095269,2.101756,2.101928, 2.101952],
[2.599545,3.063189,3.081057,3.083965,3.085387]],[[1.369579,1.462782,1.465478,1.465547,1.465578],[1.920444,2.105089,2.111554,2.111633,2.111780],[2.622261,3.082044,3.099112,3.102763,3.103441]],
[[1.376458,1.469854,1.472349,1.472453,1.472531],[1.933259,2.114999,2.121438,2.121521,2.121702],[2.643282, 3.101187,3.118645,3.121287,3.121690]],[[1.384282,1.476966,1.480404,1.480559,1.480576],
[1.948464,2.126446,2.133841,2.133936,2.134194],[2.667571,3.121373,3.143574,3.145042,3.145468]],[[1.391945,1.485232,1.487427,1.487522,1.487593], [1.961357,2.138602,2.144182,2.144204,2.144253],
[2.689127, 3.146294,3.161756,3.163726,3.164096]],[[1.401008,1.492834,1.494943,1.495068,1.495171],[1.978307,2.151764,2.156617,2.156934,2.157615],[2.716153,3.169754,3.188961, 3.192023,3.193316]],
[[1.411670,1.500951,1.503252,1.503442,1.503842],[1.993349,2.165997,2.173272,2.173545,2.173771],[2.744999,3.198712,3.214096,3.216680,3.217122]],[[1.420877,1.510675,1.511970,1.511986,1.512084],
[2.010536,2.184522,2.190798,2.191140,2.191611],[2.769564,3.225578,3.237936,3.240050,3.240793]],[[1.433263,1.519837,1.521600,1.521629,1.521645],[2.031463, 2.201170,2.208535,2.208754,2.209073],
[2.799616,3.252830,3.274006,3.274860,3.276932]],[[1.443600,1.532697,1.533838,1.534082,1.534365],[2.052896,2.218109,2.224369,2.224911,2.225384],[2.842197,3.294485, 3.309482,3.310631,3.311442]],
[[1.455112,1.544403,1.545380,1.545547,1.545562],[2.073102,2.241080,2.246534,2.246710,2.247286],[2.882628,3.328245,3.339912,3.340922,3.341217]],[[1.471607,1.559106,1.560338,1.560338,1.560361],
[2.096964,2.264123,2.269144,2.269487,2.269782],[2.923850,3.368009,3.382317,3.383336,3.384313]],[[1.488860,1.575721,1.576582,1.576582,1.576732],[2.121520,2.290532,2.294708,2.295162,2.295703],
[2.971624,3.413466,3.423838,3.424968,3.425139]],[[1.507300,1.596956,1.597971,1.597971,1.597971],[2.151915,2.320520,2.325255,2.325522,2.325522],[3.029458,3.460251,3.473393,3.474227,3.474227]],
[[1.531756,1.618122,1.618397,1.618397,1.618397],[2.200397,2.355904,2.358182,2.359174,2.359174],[3.087042,3.515429,3.529342,3.529363,3.529363]],[[1.560381,1.649405,1.649405,1.649405,1.649405],
[2.247114,2.408982,2.411867,2.411867,2.411867],[3.172736,3.606243, 3.619386, 3.620481,3.620959]],[[1.604588,1.685943,1.685943,1.685943,1.685943],[2.308004,2.464537,2.465797,2.465797,2.465797],
[3.289425,3.718164,3.734348,3.736920,3.736979]],[[1.673977,1.745509, 1.745509,1.745509,1.745509],[2.434576,2.568862,2.570255,2.570255,2.570255],[3.454727,3.935357,3.941029,3.941029,3.941029]]]// 1. GET  CRITICAL VALUE FROM THE CRITICAL VALUE TABLE.
var tblindex = criticalVTable.slice(alphaIndex);
var tblind = tblindex[0];
var htablindex = tblind[h];
var criticalValue = htablindex[period];///
//MONITORING  PERIOD
//// Filter to the area of interest, mask clouds, add variables.
var MonitorLandsat = moniCollection//.map(addNDVI).map(addDependents).map(addHarmonics(harmonicFrequencies));// Do the prediction.
var predictedValues = MonitorLandsat.map(function(image) {return image.addBands(image.select(independents).multiply(harmonicTrendCoefficients).reduce('sum').rename('predicted'))
});
var mcomputResiduals = function(image){return image.addBands(image.select(dependent).subtract(image.select('predicted')).reduce('sum').rename('residual'))
}
// Calculate the MOSUM for monitoring period///...calculate the Cumsum of residuals and do a shift trick to get MOSUM
var mresiduals = predictedValues.map(mcomputResiduals);
//add 5 band
//print(mresiduals)
var mresid = mresiduals.select('residual')var moResiduals = mresid.toArray();
var histResidual =  res.toArray();
var moMoresid = ee.ImageCollection([histResidual,moResiduals]);
var MonResid = moMoresid.toArrayPerBand(imageAxis);
var moRes = MonResid.toArray()var cumSumRes = moRes.arrayAccum(imageAxis, 'sum');
var imZer0 = ee.Image(0).int();
var ksizeplu = ksize.add(1)
var dubImZer0 = imZer0.arrayRepeat(imageAxis,ksize.int());
var fx10 = dubImZer0.arrayRepeat(bandAxis, 1);
var verc= fx10.toArray().float();
var shiftedCumR = ee.ImageCollection([verc, cumSumRes]);
var ShfCumSumRight1 = shiftedCumR.toArrayPerBand(imageAxis);
var ShfCumSumRight = ShfCumSumRight1.toArray()
var shiftedCumL = ee.ImageCollection([cumSumRes, verc]);
var ShfCumSumLeft = shiftedCumL.toArrayPerBand(imageAxis);
var ShiftDifference = ShfCumSumLeft.subtract(ShfCumSumRight);
var msk = ShiftDifference.arrayMask(ShfCumSumLeft);
var ter = ee.Image(0).float().toArray();
var tesser = ter.arrayRepeat(bandAxis, 1);
var tessermaske = ee.ImageCollection([tesser, msk]);
var tessermsk = tessermaske.toArrayPerBand(imageAxis);///... now slash off unwanted first part of Cumsum of residuals to get MOSUM residuals
var cotx2 = tessermsk.arrayLength(imageAxis);
var cumKsizeDiff = cotx2.subtract(ksize).int();var knsiz1 = imZer0.arrayRepeat(imageAxis, ksize.int());
var knsiz2 = knsiz1.arrayRepeat(bandAxis, 1);
var knsiz =knsiz2.toArray();var dcumn = ee.Image(1).multiply(cumKsizeDiff);
var mask =dcumn.select("constant").gt(0);
var dcumnMask= dcumn.updateMask(mask);
var imOnes = ee.Image(1).int();
var knm1 = imOnes.arrayRepeat(imageAxis, dcumnMask);//bugAndDebug
var knm2 = knm1.arrayRepeat(bandAxis, 1);
var knm = knm2.toArray();var shmasker = ee.ImageCollection([knsiz, knm]);var Shfmask = shmasker.toArrayPerBand(imageAxis);
var moResmosum = tessermsk.arrayMask(Shfmask);var dumoSigmaStanda = sigmaStanda.arrayRepeat(bandAxis,1);
var mobleng = moResmosum.arrayLength(imageAxis);var moSigmaStanda = dumoSigmaStanda.arrayRepeat(imageAxis,mobleng).toArray();///... now standardise the mosum
var moResmosumStanda1 = moResmosum.divide(moSigmaStanda);
var moResmosumStanda = moResmosumStanda1.abs();
///...Time stamp...needed for time of change later// Number of milliseconds in a day.
var MSEC_PER_DAY = 365.25 * 24 * 3600 * 1000;function makeDate(image) {// Create a band containing the time of each image,// converted to days since the epoch (1/1/1970).var md = image.metadata("system:time_start").divide(MSEC_PER_DAY).add(1970);// Copy the mask from band 0.var mask = image.select(0).mask();return md.mask(mask);
}var vsr = MonitorLandsat.map(makeDate);var time = vsr.toArray();///...magnitude ..needed for magnitude change latervar magnTudmas = moResiduals;/// 3.Calculate the critical borders///.. generating a sequence of values  from N obs in historical sample to N obs in monitoring sample with //interval of 1 interval and this is not a straight forward procedurevar monisize1 = moResiduals.arrayLength(imageAxis);
var signb = ee.Image(1).toArray()
var kplushistzie = signb.add(histtimeseiesLength);//logPlus
var fo1 =  kplushistzie.divide(histtimeseiesLength);
var foLogplus1 = fo1.log().max(1);//border
var foLop1 = foLogplus1.multiply(2);
var foLopSq1 = foLop1.sqrt();
var foLopSqx = foLopSq1.arrayRepeat(imageAxis,monisize1.int());var criticalBorder1 = foLopSqx.multiply(ee.Image(criticalValue));
var criticalBorder = criticalBorder1.arrayRepeat(bandAxis, 1);///...now create masks to mask out parts that belong to the history periodvar yeOnes = ee.Image(1).int();
//var timeMa = time.arrayLength(imageAxis)
var moRepeat1 = yeOnes.arrayRepeat(imageAxis, monisize1.int()).toArray();var moRepeat  = moRepeat1.arrayRepeat(bandAxis, 1);
var rMa = moResmosumStanda.arrayLength(imageAxis);
var rt = rMa.subtract(monisize1);
var rtx = ee.Image(0).int();
var mask =rt.select("array").gt(0);
var rtMask= rt.updateMask(mask);
var rhRepeat1 = rtx.arrayRepeat(imageAxis, rtMask.int()).toArray();//bugAndDebug2var rhRepeat  = rhRepeat1.arrayRepeat(bandAxis, 1);var rZOnesx = ee.ImageCollection([rhRepeat,moRepeat]);
var rMZOnesxo = rZOnesx.toArrayPerBand(0);
var rMZexOnesx = rMZOnesxo.arrayRepeat(bandAxis, 1);/// ...now start masking
var xmoResmosumStanda = moResmosumStanda.arrayMask(rMZexOnesx);
var  chiefMasker = xmoResmosumStanda.divide(criticalBorder).int();/// 4. Check if a break is detected///... Now mask out the positions where critical boundary is not closed
// you do this for change, time and magnitude of change image
var change = xmoResmosumStanda.arrayMask(chiefMasker);
var timechange = time.arrayMask(chiefMasker);
var magniofchange = magnTudmas.arrayMask(chiefMasker);
///... now slice out the first position the change has occured
var firstChange = change.arraySlice(0,0,1).toArray();
var timeOfchange = timechange.arraySlice(0,0,1).toArray();
var magnitude = magniofchange.arraySlice(0,0,1).toArray();
var criticalBorder1 = criticalBorder.arraySlice(0,0,1).toArray();///...mask out locations where no change is detected;
var time0x = ee.Image(0).int().toArray();//nobug
var timZx = time0x.arrayRepeat(bandAxis,1);//nobugvar mas = ee.ImageCollection([magnitude,timZx]);
var timOChx = ee.ImageCollection([timeOfchange,timZx]);
var firstChx = ee.ImageCollection([firstChange,timZx]);var timeOfCvx = timOChx.toArrayPerBand(0);
var tMa = mas.toArrayPerBand(0);//bug
var firsttCha = firstChx.toArrayPerBand(0);var cox = tMa.arrayReduce('sum', [imageAxis]);
var Timecox = timeOfCvx.arrayReduce("max", [imageAxis]);
var firsttChab = firsttCha.arrayReduce('max', [imageAxis]);var Cnk1 = cox.mask(cox.arrayGet([0,0]).neq(0)).toArray();
var timeCnk1 = Timecox.mask(Timecox.arrayGet([0,0]).neq(0)).toArray();
var firtsChnk = firsttChab.mask(firsttChab.arrayGet([0,0]).neq(0)).toArray();
var firstCha = firtsChnk.subtract(criticalBorder1);///... flatten the array to allow for use of palletes
var timeCnk2 = timeCnk1.arrayFlatten([['x'],['y']]);
var Cnk = Cnk1.arrayFlatten([['x'],['y']]);
var firstChaDif = firstCha.arrayFlatten([['x'],['y']]);
timeCnk2 = timeCnk2.select(['.*'],['breakTime']).set({'bfast:label' : 'Time of change', 'bfast:result' : 'breakTime'})
Cnk = Cnk.select(['.*'],['breakMagnitude']).set({'bfast:label' : 'Magnitude of change', 'bfast:result' : 'breakMagnitude'})//4.导出图像，名字最好根据数据集和图像类型来，图像类型分为改变的程度和改变的时间
Export.image.toDrive({image:timeCnk2,//turn the namedescription: 'timeSnowCoverECMWF_ERA5_LAND_MONTHLY',//turn the name tooscale: 100,//region:roi,//maxPixels:1e13,fileFormat: 'GeoTIFF',});
Map.addLayer(timeCnk2, {min:yearStart, max:yearEnd,'palette' : '00BFFF,CC2EFA,A901DB,6A0888,5858FA,0101DF,2E2EFE,0B0B61'},"Time of change");
Map.addLayer(Cnk,{palette: ['blue', 'green', 'red']}, 'Magnitude of change');
print("done")

博主选择"ECMWF/ERA5_LAND/MONTHLY"数据集的snow_cover雪覆盖数据集进行计算，最后得到的改变程度和改变时间分别如下所示:

在Arcgis里可视化的变化程度和变化时间分别为：

博主也是刚开始用GEE，所以这个改编仍然有不足之处：在改编的时候，总是出现Tile error: Copies cannot be negative的报错，在历经磨难的排查bug后，最后发现是arrayrepeat这个函数不接受负值，原算法和博主导入的NDWI数据在计算时没有问题，而换成GEE上的数据集后很多都会报这样的错,最后博主选择掩膜掉负值，然后发现被掩码的地方恰好是湖泊这种没有我选择的雪覆盖的地方，于是我猜测出现负值的原因很可能是原数据集缺失的地方，不过最终结果准确度尚可。

基于GEE的bfastmonitor的改编相关推荐

基于GEE与哨兵1号影像数据提取水体
在进行光学数据提取水体时,常会发现部分时间.部分区域内由于云的存在而出现大面积的空白区域,从而使得我们提取水体面积过程中存在不精确的问题.因此,基于微波数据进行云的提取就成为我们较好的一个选择,通过阅 ...
【基于GEE计算年度植被覆盖度】
基于GEE计算年度植被覆盖度 1.植被覆盖度计算公式参考地理学报的<2001-2013年华北地区植被覆盖度与干旱条件的相关分析> 2.代码 var roi=ee.FeatureCollec ...
基于GEE平台提取水体
本文主要介绍如何利用GEE平台与哨兵2号影像提取水体.水体的提取主要是基于NDWI指数进行,当然,通过改变波段运算,也可以根据MNDWI进行提取.主要代码如下: 下面展示一些内联代码片. //首先对 ...
基于GEE平台土地类型分类
GEE可以进行大尺度的地物类型分类,主要原理就是根据不同地物类型的光谱反射率的不同,目前地物分类的算法有很多种,比如随机森林算法等等,本文章主要基于随机森林算法进行大尺度地物分类,代码如下所示: va ...
基于GEE平台的植被覆盖度（FVC）像元二分法计算
一.植被覆盖度计算方法植被覆盖度FVC(Fractional Vegetation Cover)定义为单位面积内绿色植被冠层垂直投影面积所占比例.FVC是衡量地表植被状况的重要指标之一,也是区域生态 ...
GEE学习记录（一）基于GEE利用LANDSAT 8数据计算遥感生态指数（RSEI）
最近老师让看一下关于GEE的东西,实现大面积的反演.计算地表温度等,也算熟悉一下.参考网上很多大佬的文章,按照自己的思路和想法算出了RSEI,参考的文章都有列出来. 目录所用数据集影像数据矢量数 ...
GEE：基于GEE的单个湖泊的实时水体提取（以武汉东湖为例）
前言博主不主修遥感方向,不是专业人士,由于毕设需要使用GEE,故临时学习并做了记录,由于是博主自己钻研的,不知道有无其他更便捷的方式.若文中有错误欢迎指正. 在做毕设的时候是对武汉市所有湖泊进行分别 ...
基于GEE黑龙江省大宗农作物空间分布（注释+全套代码）
描述:利用gee提取黑龙江省主要作物大豆,玉米,水稻,小麦,2001.2006.2011.2016.2021年这五年的黑龙江省农作物空间分布 //第一步:设置研究区 // ############## ...
【Earth Engine】基于GEE对季节性地物进行分类（多源数据叠图+监督分类）
目录 1 简介与摘要 2 采样点的选取 3 合成多季节多波段影像 4 分类器参数设置 5 监督分类的实现与分类精度计算 6 影像的显示与结果 7 后记 1 简介与摘要最近导师给了个新方向,其中要用到 ...
基于GEE的制作全球任意地方时间序列数据动画的方法
大家好,我是南南今天来教大家玩个好东西(超简单) 众所周知,由于卫星遥感观测具有重访性特点,迄今已经积累了大量的各种地表参数遥感时间序列产品,这些时间序列数据较为真实地反映了地表在一个长时间范围内的 ...

基于GEE的bfastmonitor的改编

基于GEE的bfastmonitor的改编相关推荐

最新文章

热门文章