1.概略

本記事では、GEE (Javascript)でCCDCを実行するコードを紹介します。
GEEにはCCDCのAPIが実装されていますが、なぜかエラーが出力されることが多く、使えなかったため自分で組んでみます。必要なところだけでAPIを使用します。
コードを組みにあたり、以下のサイトを参考にしました。
・CCDC GUI Implementation
・CCDC github
・CCDC API Reference
・CCDC API Ducumentation

本記事で紹介するコードのワークフローは、
1.CCDCのChange Detectionの出力
2.CCDCのClassificationの学習データの作成
3.CCDCのClassificationの出力
となっています。
全ての工程を一つのプログラムにまとめるとmemory不足でエラーが出力されるため分けています。

また、本記事ではCCDCの内容について詳細に説明していません。

2.コードの実装

2.0 共通のコード

CCDCは時系列Landsat画像を高密度で処理しています。APIを使用すると簡単に入力データを得ることができますが、研究などで画像に様々な加工を施したり、新機体の導入やデータの更新に対応できないなどのデメリットがあります。
そこで、APIに頼らない高密度時系列Landsatの取得するためのコードを紹介します。
(CCDC githubの一部を変更して、そのまま持ってきたものですが)
異なる機体のLandsatを用いるため、DN値を統一するためにも不可欠なものです。
APIを用いた際に生じたエラーの原因がRemove Cloudなどのピクセルの選別の部分だったため、変更しています。また、TEMPというバンドを除去しています。
以下は関数です。

// Mask for Landsat 4, 5
function prepareL4L5(image){
  var bandList = ['SR_B1', 'SR_B2','SR_B3','SR_B4','SR_B5','SR_B7']
  var nameList = ['BLUE', 'GREEN', 'RED', 'NIR', 'SWIR1', 'SWIR2']
  var scaling = [10000, 10000, 10000, 10000, 10000, 10000]
  var scaled = ee.Image(image).select(bandList).rename(nameList).divide(ee.Image.constant(scaling))
  
  var qa = image.select("QA_PIXEL");
  var maskClouds = qa.bitwiseAnd(1 << 3).eq(0); // Remove Cloud
  var maskShadows = qa.bitwiseAnd(1 << 4).eq(0); // Remove Shadow
  var maskSnows = qa.bitwiseAnd(1 << 5).eq(0); // Remove Snow
  // Gat valid data mask, for pixels without band saturation
  var mask2 = image.select('QA_RADSAT').eq(0);
  var mask3 = image.select(bandList).reduce(ee.Reducer.min()).gt(0);
  // Mask hazy pixels. Aggressively filters too many images in arid regions (e.g Egypt)
  // unless we force include 'nodata' values by unmasking
  var mask4 = image.select("SR_ATMOS_OPACITY").unmask().lt(300);
  return ee.Image(image).addBands(scaled).updateMask(maskClouds.and(maskShadows).and(maskSnows).and(mask2).and(mask3).and(mask4));
}
// Mask for Landsat 7
function prepareL7(image){
  var bandList = ['SR_B1', 'SR_B2','SR_B3','SR_B4','SR_B5','SR_B7']
  var nameList = ['BLUE', 'GREEN', 'RED', 'NIR', 'SWIR1', 'SWIR2']
  var scaling = [10000, 10000, 10000, 10000, 10000, 10000]
  var scaled = ee.Image(image).select(bandList).rename(nameList).divide(ee.Image.constant(scaling))
  
  var qa = image.select("QA_PIXEL");
  var maskClouds = qa.bitwiseAnd(1 << 3).eq(0); // Remove Cloud
  var maskShadows = qa.bitwiseAnd(1 << 4).eq(0); // Remove Shadow
  var maskSnows = qa.bitwiseAnd(1 << 5).eq(0); // Remove Snow
  // Gat valid data mask, for pixels without band saturation
  var mask2 = image.select('QA_RADSAT').eq(0)
  var mask3 = image.select(bandList).reduce(ee.Reducer.min()).gt(0)
  // Mask hazy pixels. Aggressively filters too many images in arid regions (e.g Egypt)
  // unless we force include 'nodata' values by unmasking
  var mask4 = image.select("SR_ATMOS_OPACITY").unmask().lt(300)
  // Slightly erode bands to get rid of artifacts due to scan lines
  var mask5 = ee.Image(image).mask().reduce(ee.Reducer.min()).focal_min(2.5)
  return ee.Image(image).addBands(scaled).updateMask(maskClouds.and(maskShadows).and(maskSnows).and(mask2).and(mask3).and(mask4).and(mask5))
}
// Mask for Landsat 8, 9
function prepareL8L9(image){
  var bandList = ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7']
  var nameList = ['BLUE', 'GREEN', 'RED', 'NIR', 'SWIR1', 'SWIR2']
  var scaling = [10000, 10000, 10000, 10000, 10000, 10000]

  var validTOA = [66, 68, 72, 80, 96, 100, 130, 132, 136, 144, 160, 164]
  var validQA = [322, 386, 324, 388, 836, 900]

  var scaled = ee.Image(image).select(bandList).rename(nameList).divide(ee.Image.constant(scaling))
  var qa = image.select("QA_PIXEL");
  var maskClouds = qa.bitwiseAnd(1 << 3).eq(0); // Remove Cloud
  var maskShadows = qa.bitwiseAnd(1 << 4).eq(0); // Remove Shadow
  var maskSnows = qa.bitwiseAnd(1 << 5).eq(0); // Remove Snow
  var mask2 = image.select('QA_RADSAT').eq(0)
  var mask3 = image.select(bandList).reduce(ee.Reducer.min()).gt(0)
  var mask4 = ee.Image(image).select(['SR_QA_AEROSOL']).remap(validTOA, ee.List.repeat(1, validTOA.length), 0)
  return ee.Image(image).addBands(scaled).updateMask(maskClouds.and(maskShadows).and(maskSnows).and(mask2).and(mask3).and(mask4))
}

// Define Function for Vegetable Index
function calcNDVI(img){
  return ee.Image(img).normalizedDifference(["NIR", "RED"]).rename("NDVI")
}
function calcNBR(image) {
  return ee.Image(image).normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
function tcTrans(image) {

    // Calculate tasseled cap transformation
    var brightness = image.expression(
        '(L1 * B1) + (L2 * B2) + (L3 * B3) + (L4 * B4) + (L5 * B5) + (L6 * B6)',
        {
            'L1': image.select('BLUE'),
            'B1': 0.2043,
            'L2': image.select('GREEN'),
            'B2': 0.4158,
            'L3': image.select('RED'),
            'B3': 0.5524,
            'L4': image.select('NIR'),
            'B4': 0.5741,
            'L5': image.select('SWIR1'),
            'B5': 0.3124,
            'L6': image.select('SWIR2'),
            'B6': 0.2303
        });
    var greenness = image.expression(
        '(L1 * B1) + (L2 * B2) + (L3 * B3) + (L4 * B4) + (L5 * B5) + (L6 * B6)',
        {
            'L1': image.select('BLUE'),
            'B1': -0.1603,
            'L2': image.select('GREEN'),
            'B2': -0.2819,
            'L3': image.select('RED'),
            'B3': -0.4934,
            'L4': image.select('NIR'),
            'B4': 0.7940,
            'L5': image.select('SWIR1'),
            'B5': -0.0002,
            'L6': image.select('SWIR2'),
            'B6': -0.1446
        });
    var wetness = image.expression(
        '(L1 * B1) + (L2 * B2) + (L3 * B3) + (L4 * B4) + (L5 * B5) + (L6 * B6)',
        {
            'L1': image.select('BLUE'),
            'B1': 0.0315,
            'L2': image.select('GREEN'),
            'B2': 0.2021,
            'L3': image.select('RED'),
            'B3': 0.3102,
            'L4': image.select('NIR'),
            'B4': 0.1594,
            'L5': image.select('SWIR1'),
            'B5': -0.6806,
            'L6': image.select('SWIR2'),
            'B6': -0.6109
        });

    var bright =  ee.Image(brightness).rename('BRIGHTNESS');
    var green = ee.Image(greenness).rename('GREENNESS');
    var wet = ee.Image(wetness).rename('WETNESS');
    var tasseledCap = ee.Image([bright, green, wet])
    return tasseledCap
}
function calcIndices(collection) {
  return collection.map(function(image) {
    var NDVI =  calcNDVI(image)
    var NBR = calcNBR(image)
    var TC = tcTrans(image)
    return image.addBands([NDVI, NBR, TC])
  })
}
var getVal = function(panel, index) {
  return panel.widgets().get(index).getValue()}

var checkParams = function() {
    var crit1 = PROPS.region.geoType != ''
    var crit2 = PROPS.training.trainingStrategy != ''
    var crit3 = PROPS.classifier.classifier
    var crit4 = PROPS.predictorsLoaded
    if (crit1 && crit2 && crit3 && crit4) {
      WIDGS.main.runButton
        .widgets()
        .get(0)
        .setDisabled(false)
  }
}

以下は事前に定義する引数です。

var point = ee.Geometry.Point(lon, lat); //任意の対象地
var area = point.buffer(40000).bounds()
Map.centerObject(area, 7);
Map.addLayer(area);
// Define Study Period
var start = "2000-01-01";
var end = "2024-12-31";
// Define Bands
var bands = ["GREEN", "NDVI", "WETNESS", "SWIR1"];
var segments = 10;
var segs = ["S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9", "S10"];
var detectionBands = ["NDVI", "WETNESS"]
// want date
var wantDate = "2000-12-01";
// CCDC API
var utils = require('users/parevalo_bu/gee-ccdc-tools:ccdcUtilities/api');
// 使用するCoefficients
var SELECT_COEFS = ["INTP", "SLP", "COS", "SIN", "COS2", "SIN2", "COS3", "SIN3", "RMSE"]
// Get ImageCollection
var collection5 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').filterDate(start, end).map(prepareL4L5);
var collection7 = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2').filterDate(start, end).map(prepareL7);
var collection8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2').filterDate(start, end).map(prepareL8L9);
var collection9 = ee.ImageCollection('LANDSAT/LC09/C02/T1_L2').filterDate(start, end).map(prepareL8L9);
var col = collection5.merge(collection7).merge(collection8).merge(collection9);
// Define Date Information
var dateParams = {inputFormat: 3, inputDate: wantDate, outputFormat: 1}
var formattedDate = utils.Dates.convertDate(dateParams)

2.1 CCDCのChange Detectionの出力

ここでは、対象エリアの設定、CCDCのChange Detectionの部分の画像を出力するまでのコードを紹介します。

col = col.filterBounds(area);
var indices = calcIndices(col).select(bands);
// CCDC
var result = ee.Algorithms.TemporalSegmentation.Ccdc({
  collection: indices, 
  breakpointBands: ["NDVI", "WETNESS", "GREEN", "SWIR1"],
  tmaskBands: ["GREEN", "SWIR1"],
  lambda: 20/10000,
  minObservations: 5
})
// All data
var converted = utils.CCDC.buildCcdImage(result, segs.length, detectionBands);
// Multi Coefficients
var coefs = utils.CCDC.getMultiCoefs(converted, formattedDate, detectionBands, SELECT_COEFS, true, segs, 'after')
// Obtain synthetic image
var synt = utils.CCDC.getMultiSynthetic(converted, formattedDate, 1, detectionBands, segs)
// Obtain Change Information
var startParams = {inputFormat: 3, inputDate: start, outputFormat: 1}
var endParams = {inputFormat: 3, inputDate: end, outputFormat: 1}
var formattedStart = utils.Dates.convertDate(startParams)
var formattedEnd = utils.Dates.convertDate(endParams)
var mag = utils.CCDC.filterMag(converted, formattedStart, formattedEnd, detectionBands[0], segs);
mag = mag.float()
// Export
Export.image.toDrive({image:converted, description:"CCDC_Result", region:area, scale:30, crs:"EPSG:4326"});
Export.image.toDrive({image:coefs, description:"CCDC_Coefficients_"+wantDate, region:area, scale:30, crs:"EPSG:4326"});

2.2 CCDCのClassificationの学習データの作成

次に、Classificationを実行するために必要な学習データを作成するコードを紹介します。
素となる学習データは、APIとgee-ccdcのサイトに載っているデータが使えなかったためusers/bullocke/ccdc/public/EastAfrica_Training_Example_Dataを使用しました。
CCDC GUI Implementationで使用されていたものです。

// 対象地とTraining Dataの地域は異なる。
// 対象地は任意、Training Dataはアフリカ
var trainingPathString = 'users/bullocke/ccdc/public/EastAfrica_Training_Example_Data';
var trainingData = ee.FeatureCollection(trainingPathString);
trainingData = trainingData.select(["landcover", "year"]);
var polygon = trainingData.bounds();
// ImageCollection
col = col.filterBounds(polygon);
var indicesForClassification = calcIndices(col).select(bands) ;
// CCDC
var resultForClassification = ee.Algorithms.TemporalSegmentation.Ccdc({
  collection: indicesForClassification,  breakpointBands: ["NDVI", "WETNESS", "GREEN", "SWIR1"], 
  tmaskBands: ["GREEN", "SWIR1"], lambda: 20/10000, minObservations: 5
})
var ccdImage = utils.CCDC.buildCcdImage(resultForClassification, segs.length, detectionBands).clip(polygon);
// Additional Data
var demImage = ee.Image('USGS/SRTMGL1_003').rename('ELEVATION');
var slope = ee.Terrain.slope(demImage).rename('DEM_SLOPE');
var aspect = ee.Terrain.aspect(demImage).rename('ASPECT');
var bio = ee.Image('WORLDCLIM/V1/BIO').select(['bio01','bio12']).rename(['TEMPERATURE','RAINFALL']);
var ancillary = ee.Image.cat([demImage, slope, aspect, bio]);

// Make Dataset
var coefForTraining = utils.CCDC.getMultiCoefs(ccdImage, formattedDate, detectionBands, SELECT_COEFS, true, segs,'after').addBands(ancillary).unmask();
coefForTraining = coefForTraining.clip(polygon);
var withCoefs = coefForTraining.sampleRegions({collection: trainingData, scale: 30, tileScale: 16, geometries: true});
var bandsToClassify = ee.List(detectionBands.map(function(b) {return SELECT_COEFS.map(function(c) {return b.concat('_').concat(c)})})).flatten();
var savedCoefs = trainingData.first().propertyNames().containsAll(bandsToClassify);
withCoefs = ee.FeatureCollection(ee.Algorithms.If(savedCoefs, trainingData, withCoefs));
// Export
Export.table.toDrive({collection: withCoefs, description:'CCDC_TrainingData',fileFormat: 'CSV'});

Google Driveに出力されたcsvファイルをダウンロードします。これはQGISでも使用することができます。
ただし、地物として出力する際のフォーマットをshpファイルとすると、格納できる列名の長さが限られいるためgeojsonがおすすめです。
出力されたcsvファイルをassetsからGEEの再入力します。assetsのNewを押すと以下が出てきます。

これのCSV file(.csv)を押して、ダウンロードしたcsvファイルを選択します。
これで、2.3の準備ができました。

2.3 CCDCのClassificationの出力

コードを実行する前に、インポートしたcsvファイルをコードに入力します。Assetsにインポートしたファイルの→を押すとコード内にデータを取り込むことができます。
取り込んだデータのTable1クリックするとファイル名を変更することができます。
今回のファイル名はTrainingDataとします。

col = col.filterBounds(area);
var indices = calcIndices(col).select(bands) ;
// CCDC
var ccdcResult = ee.Algorithms.TemporalSegmentation.Ccdc({
  collection: indices,  breakpointBands: detectionBands, 
  tmaskBands: ["GREEN", "SWIR1"], lambda: 20/10000, minObservations: 5
})
var ccdImage = utils.CCDC.buildCcdImage(ccdcResult, segs.length, detectionBands).clip(area);
// Additional Data
var demImage = ee.Image('USGS/SRTMGL1_003').rename('ELEVATION');
var slope = ee.Terrain.slope(demImage).rename('DEM_SLOPE');
var aspect = ee.Terrain.aspect(demImage).rename('ASPECT');
var bio = ee.Image('WORLDCLIM/V1/BIO').select(['bio01','bio12']).rename(['TEMPERATURE','RAINFALL']);
var ancillary = ee.Image.cat([demImage, slope, aspect, bio]);
// Classification
var ancillaryFeatures = ["ELEVATION","ASPECT","DEM_SLOPE","RAINFALL","TEMPERATURE"];
var coef = ["INTP", "SLP", "RMSE"];
var classifier = ee.Classifier.smileRandomForest({numberOfTrees: 150, variablesPerSplit: null,
                                                  minLeafPopulation: 1, bagFraction: 0.5, maxNodes: null});                                 
var classificationResult = utils.Classification.classifySegments(
    ccdImage, segs.length, detectionBands, ancillary, ancillaryFeatures,
    trainingData, classifier, null, "landcover", coef).clip(area).toUint8();
print(classificationResult);
Map.addLayer(classificationResult.select(0).randomVisualizer(), {}, 'Seg1 Classification')
Export.image.toDrive({image:classificationResult, description: "ccdcClassification", scale:30, crs:"EPSG:4326"});

九州の霧島を対象としたClassificationの結果が以下です。

遠目に見るとよさげに見えますが、近づくとゴマ塩ノイズ (塊となっていない1ピクセル)が結構あるため、後処理は不可欠です。

3.終わりに

今回は、CCDCをAPIを利用しつつもデータの取得や出力を自作のコードで実装しました。
APIを利用するとエラーが出ていましたが、それらを解消することができました。
しかし、ゴマ塩ノイズの除去や精度検証、およびCCDCの説明は行っていません。不完全なコードではありますが、誰かの参考になればいいなと思います (ノイズの除去はもしかしたら、記事を作成するかもしれません)。