《gee大神之路》月下醉吟

13、十三 第1节：通

　　第1节：通过 CCDC 了解时间分割

　　通过测试反射率时间序列中的结构断裂来检测像素级的光谱变化。在 Earth Engine 中，此过程称为“时间分割”，因为像素级时间序列根据唯一反射率周期进行分割。它通过将调和回归模型拟合到时间序列中的所有光谱带来实现这一点。模型拟合从时间序列的开头开始，并以“在线”方法及时向前推进以进行变化检测。这些系数用于预测未来的观测值，如果未来观测值的残差超过多次连续观测值的统计阈值，则算法会标记已发生变化。更改后，将拟合新的回归模型，并且该过程将持续到时间序列结束。

　　加载 CCDC 界面你将能够导航到任何位置，选择要绘制的 Landsat 光谱带或索引，然后单击地图以查看 CCDC 在你所在位置的拟合情况。点击。在本次练习中，我们将研究巴西朗多尼亚州的景观动态。我们可以使用左下角的面板输入以下坐标（纬度，经度）：-9.0002，-62.7223。将在该位置添加一个点，并且地图将放大到该点。到达那里后，单击该点并等待底部的图表加载。

　　此示例显示第一个短波红外 (S WIR1)的 Landsat 时间序列波段（如蓝点）和时间段（如彩色线）使用 CCDC 默认参数运行。第一部分代表稳定的森林，在2006年中期突然被砍伐。该算法检测到这一变化事件并随后拟合一个新的片段，代表农业的新时间模式。随着新片段的安装，检测到其他后续模式，这些新片段可能对应于收获和再生周期或不同的作物。要研究随时间变化的动态，你可以单击图表中的点，它们对应的 Landsat 图像将根据为左侧面板中的 RGB 组合选择的可视化参数添加到地图中。目前，该面板中所做的更改不会立即生效，而是必须在单击地图之前进行设置。

　　特别注意每个细分市场的特点。例如，查看每个段的平均表面反射率值。明显斜坡的存在可能表明植被再生或退化等现象。每个片段中使用的谐波数量可以代表植被的季节性（自然的或由于农业实践）或景观动态（例如，季节性洪水）。

　　第二节运行

　　上面显示的工具对于理解特定点的时间动态非常有用。然而，我们可以通过首先对一组像素运行 CCDC 算法来对更大的区域进行类似的分析。Earth Engine 中的 CCDC 函数可以采用任何ImageCollection ，最好是噪声很少或没有噪声的 ImageCollection ，例如云和云阴影已被屏蔽的Landsat ImageCollection。CCDC 包含一个内部云屏蔽算法，对于漏云具有很强的鲁棒性，但数据越干净越好。为了简化流程，我们开发了一个函数库，其中包含用于生成输入数据和处理 CCDC 结果的函数。将此行代码粘贴到新脚本中：

　　var utils=require(
　　\'users/parevalo_bu/gee-ccdc-tools:ccdcUtilities/api\');

　　对于当前练习，我们将获得 Landsat 4、5、7和8数据的 ImageCollection（集合2第1层），该数据已过滤云、云阴影、雾霾和辐射饱和像素。如果我们手动执行此操作，我们将检索每个卫星的每个ImageCollection ，应用相应的过滤器，然后将它们全部合并到一个ImageCollection中。相反，为了简化该过程，我们将使用函数getLandsat （包含在我们实用程序的“输入”模块中），然后将生成的ImageCollection过滤到2000年至2020年期间的一个小型研究区域。

　　getLandsat 函数将检索所有表面反射率波段（重命名并缩放为实际表面反射率单位）以及其他植被指数。为了简化练习，我们将仅选择要使用的表面反射带，并将以下代码添加到脚本中：

　　var studyRegion=ee.Geometry.Rectangle([
　　[-63.9533, -10.1315],
　　[-64.9118, -10.6813]
　　]);

　　// Define start, end dates and Landsat bands to use.
　　var startDate=\'2000-01-01\';
　　var endDate=\'2020-01-01\';
　　var bands=[\'BLUE\', \'GREEN\', \'RED\', \'NIR\', \'SWIR1\', \'SWIR2\'];

　　// Retrieve all clear, Landsat 4, 5, 7 and 8 observations (Collection 2, Tier 1).
　　var filteredLandsat=utils.Inputs.getLandsat({
　　collection: 2
　　})
　　.filterBounds(studyRegion)
　　.filterDate(startDate, endDate)
　　.select(bands);

　　print(filteredLandsat.first());

　　准备好ImageCollection后，我们可以指定 CCDC 参数并运行算法。对于本练习，我们将使用默认参数，这些参数在大多数情况下都工作得相当好。我们将修改的唯一参数是断点带、日期格式和 lambda。我们将在字典中设置所有参数值，并将其传递给 CCDC 函数。对于断裂检测过程，我们使用除蓝色和表面温度带（分别为\'BLUE\'和\'TEMP\' ）之外的所有带。minObservations默认值6表示标记更改所需的连续观察数。chiSquareProbability 和minNumOfYearsScaler 默认参数0.99和1.33分别控制算法检测变化的灵敏度以及检测变化所需的迭代曲线拟合过程。

　　我们将日期格式设置为 =1，它对应于小数年份，并且往往更容易解释。例如，在2010年中午检测到的变化将作为2010.5存储在像素中。最后，我们使用lambda的默认值20，但我们对其进行缩放以匹配输入的比例（表面反射率单位），并指定 maxIterations值为10000，而不是默认值25000，这可能需要更长的时间完全的。这两个参数控制曲线拟合过程。

　　为了完成输入参数，我们指定要使用的ImageCollection ，这是我们在前面的代码部分中导出的。添加以下代码：
　　// Set CCD params to use.
　　var ccdParams={
　　breakpointBands: [\'GREEN\', \'RED\', \'NIR\', \'SWIR1\', \'SWIR2\'],
　　tmaskBands: [\'GREEN\', \'SWIR2\'],
　　minObservations: 6,
　　chiSquareProbability: 0.99,
　　minNumOfYearsScaler: 1.33,
　　dateFormat: 1,
　　lambda: 0.002,
　　maxIterations: 10000,
　　collection: filteredLandsat
　　};

　　// Run CCD.
　　var ccdResults=ee.Algorithms.TemporalSegmentation.Ccdc(ccdParams);
　　print(ccdResults);

　　请注意，输出ccdResults 包含大量波段，其中一些对应于二维数组。我们将在下一节中更多地探讨这些频段。针对多个像素以交互方式运行算法的过程可能很快就会对系统造成很大的负担，从而导致内存错误。为了避免出现此类问题，我们通常首先将结果导出到 Earth Engine 资产，然后检查该资产。这种方法确保了 CCDC 成功完成运行，也让我们以后可以轻松访问结果。在本章的以下部分中，我们将使用预先计算的资产，而不是要求你自己导出资产。供你参考，导出CCDC结果所需的代码如下所示，将标志设置为 false 可以帮助你记住现在不要导出结果，而是在以下部分中使用预先计算的资源。
　　var exportResults=false
　　if (exportResults){
　　// Create a metadata dictionary with the parameters and arguments used.
　　var metadata=ccdParams;
　　metadata[\'breakpointBands\']=
　　metadata[\'breakpointBands\'].toString();
　　metadata[\'tmaskBands\']=metadata[\'tmaskBands\'].toString();
　　metadata[\'startDate\']=startDate;
　　metadata[\'endDate\']=endDate;
　　metadata[\'bands\']=bands.toString();

　　// Export results, assigning the metadata as image properties.
　　//
　　Export.image.toAsset({
　　image: ccdResults.set(metadata),
　　region: studyRegion,
　　pyramidingPolicy:{
　　\".default\": \'sample\'
　　},
　　scale: 30
　　});
　　}
　　请注意上面的元数据变量。这并不是导出每像素 CCDC 结果的严格要求，但它允许我们通过将此信息作为元数据附加到图像来记录运行的重要属性。此外，我们创建的一些与 CCDC 输出交互的工具使用用户创建的元数据来促进资产的使用。另请注意，将PyramidingPolicy的值设置为“sample”可确保输出中的所有频段都具有正确的策略。

　　作为一般规则，如果可能，请尝试使用预先存在的 CCDC 结果，如果你想尝试在本实验室练习之外自行运行它，请从非常小的区域开始。例如，本练习中的研究区域平均需要大约30分钟才能导出，但较大的图块可能需要几个小时才能完成，具体取决于集合中图像的数量和使用的参数。

　　第 3 节.提取中断信息
　　我们现在将开始探索上一节中提到的预导出的 CCDC 结果。我们将使用第三方模块调色板，简化了可视化调色板的使用。将以下代码粘贴到新脚本中：

　　var palettes=require(\'users/gena/packages:palettes\');

　　var resultsPath=
　　\'projects/gee-book/assets/F4-7/Rondonia_example_small\';
　　var ccdResults=ee.Image(resultsPath);
　　Map.centerObject(ccdResults, 10);
　　print(ccdResults);

　　第一行调用一个库，该库将有助于可视化图像。第二行包含上一节中显示的 CCDC 运行的预先计算结果的路径。
　　打印的资源将包含以下区域：

　　?tStart：每个时间段的开始日期
　　?tEnd ：每个时间段的结束日期
　　?tBreak ：检测到更改时的时间段中断日期
　　?numObs ：每个时间段使用的观察数
　　?changeProb ：表示用于变化检测的每个频段的变化概率的数值
　　?*_coefs ：输入图像集合中每个波段的回归系数
　　?*_rmse ：每个时间段和输入频段的模型均方根误差
　　?*_magnitude ：对于检测到变化的时间段，这表示变化期间的归一化残差

　　请注意，带名称和带类型旁边还有维度数（即1维、2维）。

　　这表明我们正在处理数组图像，它通常需要一组特定的函数来进行正确的操作，其中一些我们将在接下来的步骤中使用。我们将首先查看变化带，这是 CCDC 算法的关键输出之一。我们将选择包含休息时间信息的波段，并查找给定时间范围内的休息次数。

　　在同一脚本中，粘贴以下代码：
　　// Select time of break and change probability array images.
　　var change=ccdResults.select(\'tBreak\');
　　var changeProb=ccdResults.select(\'changeProb\');

　　// Set the time range we want to use and get as mask of
　　// places that meet the condition.
　　var start=2000;
　　var end=2021;
　　var mask=change.gt(start).and(change.lte(end)).and(changeProb.eq(
　　1));
　　Map.addLayer(changeProb,{}, \'change prob\');

　　// Obtain the number of breaks for the time range.
　　var numBreaks=mask.arrayReduce(ee.Reducer.sum(), [0]);
　　Map.addLayer(numBreaks,{
　　min: 0,
　　max: 5
　　}, \'Number of breaks\');

　　使用此代码，我们定义要使用的时间范围，然后生成一个掩码，该掩码将指示图像数组中在该范围内检测到的中断的所有位置，这些位置也满足变化概率为1的条件，有效地消除了一些虚假的中断。对于每个像素，我们可以计算掩模检索到有效结果的次数，表示 CCDC 检测到的中断次数。在加载图层中，看起来更亮的地方将显示更多的中断次数，可能表明从森林到农业的转变，随后是多个农业循环。请记住，检测到断裂并不总是意味着土地覆盖发生变化。自然事件、小规模扰动和季节性周期等都可能导致 CCDC 检测到中断。类似地，像素中土地覆盖状况的变化也可以被 CCDC 检测为中断，并且由于噪声时间序列或其他因素也可能发生一些错误的中断。
　　对于变化较多的地方，可视化第一次或最后一次记录中断的时间有助于了解景观中发生的变化动态。将以下代码粘贴到同一脚本中：
　　// Obtain the first change in that time period.
　　var dates=change.arrayMask(mask).arrayPad([1]);
　　var firstChange=dates
　　.arraySlice(0, 0, 1)
　　.arrayFlatten([
　　[\'firstChange\']
　　])
　　.selfMask();

　　var timeVisParams={
　　palette: palettes.colorbrewer.YlOrRd[9],
　　min: start,
　　max: end
　　};
　　Map.addLayer(firstChange, timeVisParams, \'First change\');

　　// Obtain the last change in that time period.
　　var lastChange=dates
　　.arraySlice(0, -1)
　　.arrayFlatten([
　　[\'lastChange\']
　　])
　　.selfMask();
　　Map.addLayer(lastChange, timeVisParams, \'Last change\');

　　在这里，我们使用arrayMask通过使用我们之前创建的掩码来仅保留满足条件的更改日期。我们使用函数arrayPad来填充或“填充”那些没有经历任何变化的像素，因此在tBreak band 中没有值。然后，我们选择数组中的第一个或最后一个值，并将图像从一维数组转换为常规图像，以便使用自定义调色板对其应用可视化。

　　最后，我们可以使用幅度带来可视化在我们选择的时间段内 CCDC 记录的最大变化发生的地点和时间。我们将使用 SWIR1 频段的变化幅度，以与之前相同的方式对其进行屏蔽和填充。将此代码粘贴到你的脚本中：

　　// Get masked magnitudes.
　　var magnitudes=ccdResults
　　.select(\'SWIR1_magnitude\')
　　.arrayMask(mask)
　　.arrayPad([1]);

　　// Get index of max abs magnitude of change.
　　var maxIndex=magnitudes
　　.abs()
　　.arrayArgmax()
　　.arrayFlatten([
　　[\'index\']
　　]);

　　// Select max magnitude and its timing
　　var selectedMag=magnitudes.arrayGet(maxIndex);
　　var selectedTbreak=dates.arrayGet(maxIndex).selfMask();

　　var magVisParams={
　　palette: palettes.matplotlib.viridis[7],
　　min: -0.15,
　　max: 0.15
　　};
　　Map.addLayer(selectedMag, magVisParams, \'Max mag\');
　　Map.addLayer(selectedTbreak, timeVisParams, \'Time of max mag\');

　　较深的颜色代表较新的日期，而较亮的颜色代表较旧的日期。第一个变化层显示了接近2000年时原始农业扩张的清晰模式。

　　最后一个变化层显示了同一地区最近检测到的干扰中断。图像中心的稀疏区域只有一次变化，对应于一次森林砍伐事件。未检测到中断的像素被屏蔽，因此显示下面的底图，设置为显示卫星图像。

　　我们首先取绝对值，因为幅度可以是正的也可以是负的，具体取决于变化的方向和所使用的频带。例如，SWIR1 中的正值可能表示森林损失事件，其中表面反射率从低值变为较高值。较亮的值代表该类型的事件。

　　相反，由于反射率相应下降，洪水事件将具有负值。一旦找到最大绝对值，我们就找到它在数组中的位置，然后使用该索引提取原始幅度值以及中断发生的时间。

　　第4节：手动提取系数
　　除了CCDC算法生成的变化信息之外，我们还可以将时间段的系数用于多种目的，例如土地覆盖分类。每个时间段可以被描述为具有截距、斜率和三对正弦和余弦项的调和函数，这些函数允许时间段表示在不同时间尺度上发生的季节性。这些系数以及通过比较每个预测值和实际 Landsat 值获得的均方根误差 (RMSE)是在 CCDC 算法运行时生成的。以下示例将向你展示如何检索与特定日期相交的线段的截距系数。

　　在新脚本中，粘贴以下代码：
　　var palettes=require(\'users/gena/packages:palettes\');

　　var resultsPath=
　　\'projects/gee-book/assets/F4-7/Rondonia_example_small\';
　　var ccdResults=ee.Image(resultsPath);
　　Map.centerObject(ccdResults, 10);
　　print(ccdResults);

　　// Display segment start and end times.
　　var start=ccdResults.select(\'tStart\');
　　var end=ccdResults.select(\'tEnd\');
　　Map.addLayer(start,{
　　min: 1999,
　　max: 2001
　　}, \'Segment start\');
　　Map.addLayer(end,{
　　min: 2010,
　　max: 2020
　　}, \'Segment end\');

　　检查控制台并展开打印图像信息中的条带部分。我们将使用tStart 、tEnd和SWIR1_coefs 频段，它们是数组图像，包含时间段开始的日期、日期时间段结束的日期以及 SWIR1 频段每个段的系数。运行上面的代码并将地图切换到卫星模式。使用检查器，单击图像上的任意位置，注意打印的日期数量及其多个单击像素的值。你会注意到，对于森林覆盖稳定的地方， tStart 和 tEnd通常有一个值。

　　这意味着对于那些比较稳定的地方，中央结算公司只拟合了一个时间段。另一方面，对于底图中具有可见变换的地点，日期的数量通常为两个或三个，这意味着算法分别拟合两个或三个时间段。为了简化数据的处理，我们可以选择单个段来提取其系数。粘贴以下代码并重新运行脚本：

　　// Find the segment that intersects a given date.
　　var targetDate=2005.5;
　　var selectSegment=start.lte(targetDate).and(end.gt(targetDate));
　　Map.addLayer(selectSegment,{}, \'Identified segment\');

　　在上面的代码中，我们设置了一个感兴趣的时间，在本例中是2005年中旬，然后我们找到满足在该日期之前开始和在该日期之后结束的条件的段。再次使用检查器，单击不同的位置并验证输出。满足条件的段的值为1，其他段的值0。我们可以使用此信息来选择该段的系数，使用以下代码：

　　// Get all coefs in the SWIR1 band.
　　var SWIR1Coefs=ccdResults.select(\'SWIR1_coefs\');
　　Map.addLayer(SWIR1Coefs,{}, \'SWIR1 coefs\');

　　// Select only those for the segment that we identified previously.
　　var sliceStart=selectSegment.arrayArgmax().arrayFlatten([
　　[\'index\']
　　]);
　　var sliceEnd=sliceStart.add(1);
　　var selectedCoefs=SWIR1Coefs.arraySlice(0, sliceStart, sliceEnd);
　　Map.addLayer(selectedCoefs,{}, \'Selected SWIR1 coefs\');

　　在上面的代码片段中，我们首先选择具有 SWIR1 波段系数的阵列图像。然后，使用我们之前创建的层，找到条件为真的位置，并使用它仅提取该段的系数。你可以再次使用“检查器”选项卡进行验证。
　　最后，我们现在得到的是与2005年中点相交的线段的完整系数集。系数按以下顺序排列：截距、斜率、余弦1、正弦1、余弦2、正弦2、余弦3和正弦3. 在本练习中，我们将使用以下代码提取截距系数，它是数组中的第一个元素：

　　// Retrieve only the intercept coefficient.

　　var intercept=selectedCoefs.arraySlice(1, 0, 1).arrayProject([1]);
　　var intVisParams={
　　palette: palettes.matplotlib.viridis[7],
　　min: -6,
　　max: 6
　　};
　　Map.addLayer(intercept.arrayFlatten([
　　[\'INTP\']
　　]), intVisParams, \'INTP_SWIR1\');

　　由于我们在 Landsat 表面反射率图像上运行 CCDC 算法，因此截距值应代表片段的平均反射率。但是，如果单击图像，你将看到这些值超出了0-1 范围。这是因为截距是由 CCDC 算法针对原点（例如时间0）计算的，而不是针对我们请求的年份计算的。为了检索调整后的截距以及其他系数，我们将使用不同的方法。

　　第5节：使用外部函数提取系数
　　我们在上一节中生成的代码允许我们提取单个日期的单个系数。然而，我们通常希望提取一组多个系数和频带，我们可以将其用作其他工作流程的输入，例如分类。为了简化该过程，我们将使用与我们在第2节中看到的相同的函数库。在本节中，我们将提取并可视化单个日期的不同系数，并使用同一日期的多个光谱带的截距系数生成 RGB 图像。第一步涉及确定感兴趣的日期并将 CCDC 结果从阵列图像转换为常规多波段图像，以便更轻松地操作和更快地显示。

　　在新脚本中，复制以下代码：
　　// Load the required libraries.
　　var palettes=require(\'users/gena/packages:palettes\');
　　var utils=require(
　　\'users/parevalo_bu/gee-ccdc-tools:ccdcUtilities/api\');

　　// Load the results.
　　var resultsPath=
　　\'projects/gee-book/assets/F4-7/Rondonia_example_small\';
　　var ccdResults=ee.Image(resultsPath);
　　Map.centerObject(ccdResults, 10);

　　// Convert a date into fractional years.
　　var inputDate=\'2005-09-25\';
　　var dateParams={
　　inputFormat: 3,
　　inputDate: inputDate,
　　outputFormat: 1
　　};
　　var formattedDate=utils.Dates.convertDate(dateParams);

　　// Band names originally used as inputs to the CCD algorithm.
　　var BANDS=[\'BLUE\', \'GREEN\', \'RED\', \'NIR\', \'SWIR1\', \'SWIR2\'];

　　// Names for the time segments to retrieve.
　　var SEGS=[\'S1\', \'S2\', \'S3\', \'S4\', \'S5\', \'S6\', \'S7\', \'S8\', \'S9\',
　　\'S10\'
　　];

　　// Transform CCD results into a multiband image.
　　var ccdImage=utils.CCDC.buildCcdImage(ccdResults, SEGS.length,
　　BANDS);
　　print(ccdImage);

　　在上面的代码中，我们定义了感兴趣的日期（2005-09-25）并将其转换为我们运行 CCDC 时使用的日期格式，它对应于小数年。之后，我们指定用作 CCDC 算法输入的频段。最后，我们指定分配给时间段的名称，列表长度指示每个像素检索的最大时间段数。执行此步骤是因为CCDC生成的结果存储为变长数组。例如，未检测到中断的像素将具有一个时间段，但检测到单个中断的另一像素可能具有一个或两个片段，具体取决于中断发生的时间。请求预定义的最大分段数可确保多波段图像的结构已知，并极大地方便其操作和显示。一旦我们设置了这些变量，我们调用一个函数，将结果转换为具有多个波段的图像，这些波段代表请求的段、输入波段和系数的组合。你可以在中看到图像结构控制台。

　　最后，为了提取所需频带的系数子集，我们可以使用导入库中的一个函数，称为getMultiCoefs。该函数需要以下有序参数

　　?CCDC 结果是我们在上一步中刚刚生成的多频段格式。
　　?我们想要提取系数的日期，采用 CCDC 结果运行的格式（在我们的例子中为小数年）。
　　?要检索的波段列表（即光谱波段）。
　　?要检索的系数列表，定义如下：INTP （截距）、SLP （斜率）、COS 、SIN 、COS32 、SIN2 、COS3 、SIN3和RMSE。
　　?布尔值true 或false ，指示我们是否希望针对输入日期计算截距，而不是在原点计算截距。如果为true ，则SLP 必须包含在要检索的系数列表中。
　　?段名称列表，用于在上一步中创建多波段图像。
　　?如果请求的日期没有时间段，则应用的行为：仅当日期与时间段相交时，正常才会检索值；如果没有段与日期直接相交，则 before或 after将使用紧邻请求日期之前或之后的段的值。

　　// Define bands to select.
　　var SELECT_BANDS=[\'RED\', \'GREEN\', \'BLUE\', \'NIR\'];

　　// Define coefficients to select.
　　// This list contains all possible coefficients, and the RMSE
　　var SELECT_COEFS=[\'INTP\', \'SLP\', \'RMSE\'];

　　// Obtain coefficients.
　　var coefs=utils.CCDC.getMultiCoefs(
　　ccdImage, formattedDate, SELECT_BANDS, SELECT_COEFS, true,
　　SEGS, \'after\');
　　print(coefs);

　　// Show a single coefficient.
　　var slpVisParams={
　　palette: palettes.matplotlib.viridis[7],
　　min: -0.0005,
　　max: 0.005
　　};
　　Map.addLayer(coefs.select(\'RED_SLP\'), slpVisParams,
　　\'RED SLOPE 2005-09-25\');

　　var rmseVisParams={
　　palette: palettes.matplotlib.viridis[7],
　　min: 0,
　　max: 0.1
　　};
　　Map.addLayer(coefs.select(\'NIR_RMSE\'), rmseVisParams,
　　\'NIR RMSE 2005-09-25\');

　　// Show an RGB with three coefficients.
　　var rgbVisParams={
　　bands: [\'RED_INTP\', \'GREEN_INTP\', \'BLUE_INTP\'],
　　min: 0,
　　max: 0.1
　　};
　　Map.addLayer(coefs, rgbVisParams, \'RGB 2005-09-25\');

　　对于斜率，高的正值是明亮的，而大的负值是非常暗的。剩余的森林大部分是稳定的，坡度接近于零，而经历过转变并显示出农业活动的区域往往在红色波段具有正坡度，在图像中显得明亮。类似地，对于 RMSE 图像，稳定的森林呈现出更可预测的表面反射率时间序列，这些时间序列由时间段更忠实地捕获，因此呈现出较低的 RMSE 值，在图像中显得更暗。农业地区呈现出噪声较大的时间序列，建模更具挑战性，并导致 RMSE 值更高，显得更亮。

作者有话说

　　学习使我快乐，呜呜呜呜