LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@okinasan0358(okina)

Google Earth Engine (GEE)でContinuous Change Detection and Classification (CCDC)を実装

Posted at

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を押すと以下が出てきます。
Assets_New.png
これの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の結果が以下です。
CCDC_Result.png
遠目に見るとよさげに見えますが、近づくとゴマ塩ノイズ (塊となっていない1ピクセル)が結構あるため、後処理は不可欠です。
Salt_and_Papper_Noise.png

3.終わりに

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

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?