衛星エンベディング データセットを使用した教師なし分類

GitHub で編集
問題を報告する
ページの履歴
このチュートリアルは、Satellite Embedding データセットに関するチュートリアル シリーズの一部です。概要教師あり分類回帰類似性検索もご覧ください。

前のチュートリアル(概要)では、衛星エンベディングが衛星観測と気候変数の年間の軌跡をどのようにキャプチャするかを確認しました。これにより、作物のフェノロジーをモデル化することなく、作物のマッピングにデータセットをすぐに使用できます。作物の種類をマッピングすることは、通常、作物のフェノロジーをモデル化し、その地域で栽培されているすべての作物の圃場サンプルを収集する必要があるため、難しい作業です。

このチュートリアルでは、フィールドラベルに依存せずにこの複雑なタスクを実行できる、教師なし分類アプローチによる作物マッピングを行います。この方法では、地域のローカルな知識と、世界の多くの地域で入手可能な作物の統計情報を活用します。

リージョンを選択

このチュートリアルでは、アイオワ州セロ ゴルド郡の作物タイプ地図を作成します。この郡は、米国でトウモロコシと大豆という 2 つの主要な作物が栽培されているコーンベルトにあります。このローカル知識は重要であり、モデルのキーパラメータを決定するのに役立ちます。

まず、選択した郡の境界を取得します。

// Select the region
// Cerro Gordo County, Iowa
var counties = ee.FeatureCollection('TIGER/2018/Counties');

var selected = counties
  .filter(ee.Filter.eq('GEOID', '19033'));
var geometry = selected.geometry();
Map.centerObject(geometry, 12);
Map.addLayer(geometry, {color: 'red'}, 'Selected Region', false);


図: 選択したリージョン

衛星エンベディング データセットを準備する

次に、衛星エンベディング データセットを読み込み、選択した年の画像をフィルタしてモザイクを作成します。

var embeddings = ee.ImageCollection('GOOGLE/SATELLITE_EMBEDDING/V1/ANNUAL');

var year = 2022;
var startDate = ee.Date.fromYMD(year, 1, 1);
var endDate = startDate.advance(1, 'year');

var filteredembeddings = embeddings
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.bounds(geometry));

var embeddingsImage = filteredembeddings.mosaic();

農地マスクを作成する

モデリングでは、農地以外の地域を除外する必要があります。クロップマスクの作成に使用できるグローバル データセットとリージョン データセットは多数あります。グローバルな農地データセットには、ESA WorldCover または GFSAD Global Cropland Extent Product が適しています。最近追加されたものとしては、ESA WorldCereal Active Cropland プロダクトがあります。これは、耕作地の季節ごとのマッピングです。この地域は米国にあるため、より正確な地域データセットである USDA NASS Cropland データレイヤ(CDL)を使用して、作物のマスクを取得できます。

// Use Cropland Data Layers (CDL) to obtain cultivated cropland
var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.date(startDate, endDate))
  .first();
var cropLandcover = cdl.select('cropland');
var croplandMask = cdl.select('cultivated').eq(2).rename('cropmask');

// Visualize the crop mask
var croplandMaskVis = {min: 0, max: 1, palette: ['white', 'green']};
Map.addLayer(croplandMask.clip(geometry), croplandMaskVis, 'Crop Mask');


図: 農地マスクで選択されたリージョン

トレーニング サンプルを抽出する

クロップランド マスクをエンベディング モザイクに適用します。これで、郡内の耕作地を表すピクセルがすべて残りました。

// Mask all non-cropland pixels
var clusterImage = embeddingsImage.updateMask(croplandMask);

クラスタリング モデルをトレーニングするには、衛星エンベディング画像を取得し、ランダム サンプルを取得する必要があります。対象領域にはマスクされたピクセルが多数含まれているため、単純なランダム サンプリングでは null 値のサンプルが生成される可能性があります。必要な数の null 以外のサンプルを抽出できるように、層別サンプリングを使用して、マスクされていない領域で必要な数のサンプルを取得します。

// Stratified random sampling
var training = clusterImage.addBands(croplandMask).stratifiedSample({
  numPoints: 1000,
  classBand: 'cropmask',
  region: geometry,
  scale: 10,
  tileScale: 16,
  seed: 100,
  dropNulls: true,
  geometries: true
});

サンプルをアセットにエクスポートする(省略可)

サンプルの抽出は計算コストの高いオペレーションです。抽出したトレーニング サンプルをアセットとしてエクスポートし、後続のステップでエクスポートしたアセットを使用することをおすすめします。これにより、大きなリージョンを操作する際に発生する computation timed out エラーや user memory exceeded errors を回避できます。

エクスポート タスクを開始し、完了するまで待ってから次に進みます。

// Replace this with your asset folder
// The folder must exist before exporting
var exportFolder = 'projects/spatialthoughts/assets/satellite_embedding/';

var samplesExportFc = 'cluster_training_samples';
var samplesExportFcPath = exportFolder + samplesExportFc;

Export.table.toAsset({
  collection: training,
  description: 'Cluster_Training_Samples',
  assetId: samplesExportFcPath
});

エクスポート タスクが完了したら、抽出したサンプルを特徴コレクションとしてコードに読み戻すことができます。

// Use the exported asset
var training = ee.FeatureCollection(samplesExportFcPath);

サンプルを可視化する

サンプルをインタラクティブに実行したか、フィーチャー コレクションにエクスポートしたかに関係なく、サンプル ポイントを含むトレーニング変数が作成されます。最初のサンプルを出力して検査し、トレーニング ポイントを Map に追加しましょう。

print('Extracted sample', training.first());
Map.addLayer(training, {color: 'blue'}, 'Extracted Samples', false);


図: クラスタリング用に抽出されたランダム サンプル

教師なしクラスタリングを実行する

クラスタリング モデルをトレーニングし、64D エンベディング ベクトルを任意の数の個別のクラスタにグループ化できるようになりました。Google のローカル知識に基づくと、この地域の大部分を占める 2 つの主要な作物があり、残りの部分をいくつかの他の作物が占めています。衛星エンベディングに対して教師なしクラスタリングを実行して、同様の時間的軌跡とパターンを持つピクセルのクラスタを取得できます。スペクトル特性と空間特性が類似し、フェノロジーが類似しているピクセルは、同じクラスタにグループ化されます。

ee.Clusterer.wekaCascadeKMeans() を使用すると、クラスタの最小数と最大数を指定し、トレーニング データに基づいて最適なクラスタ数を特定できます。ここで、ローカル知識がクラスタの最小数と最大数を決定するのに役立ちます。クラスタのタイプは、トウモロコシ、大豆、その他の作物のいくつかのタイプに分類されると予想されるため、クラスタの最小数を 4、クラスタの最大数を 5 に設定します。これらの数値を試して、お住まいの地域に最適な値を見つける必要がある場合があります。

// Cluster the Satellite Embedding Image
var minClusters = 4;
var maxClusters = 5;

var clusterer = ee.Clusterer.wekaCascadeKMeans({
  minClusters: minClusters, maxClusters: maxClusters}).train({
  features: training,
  inputProperties: clusterImage.bandNames()
});

var clustered = clusterImage.cluster(clusterer);
Map.addLayer(clustered.randomVisualizer().clip(geometry), {}, 'Clusters');


図: 教師なし分類で取得したクラスタ

クラスタにラベルを割り当てる

目視検査の結果、前の手順で取得したクラスタは、高解像度画像に表示されている農場の境界とほぼ一致しています。地域の知識から、2 つの最大のクラスタはトウモロコシと大豆であることがわかっています。画像の各クラスタの面積を計算してみましょう。

// Calculate Cluster Areas
// 1 Acre = 4046.86 Sq. Meters
var areaImage = ee.Image.pixelArea().divide(4046.86).addBands(clustered);

var areas = areaImage.reduceRegion({
      reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'cluster',
    }),
    geometry: geometry,
    scale: 10,
    maxPixels: 1e10
    });

var clusterAreas = ee.List(areas.get('groups'));

// Process results to extract the areas and create a FeatureCollection

var clusterAreas = clusterAreas.map(function(item) {
  var areaDict = ee.Dictionary(item);
  var clusterNumber = areaDict.getNumber('cluster').format();
  var area = areaDict.getNumber('sum')
  return ee.Feature(null, {cluster: clusterNumber, area: area})
})

var clusterAreaFc = ee.FeatureCollection(clusterAreas);
print('Cluster Areas', clusterAreaFc);

面積が最も大きい 2 つのクラスタを選択します。

var selectedFc = clusterAreaFc.sort('area', false).limit(2);
print('Top 2 Clusters by Area', selectedFc);

ただし、どのクラスタがどの作物であるかはまだわかりません。トウモロコシや大豆のフィールド サンプルがいくつかある場合は、それらをクラスタに重ねて、それぞれのラベルを特定できます。フィールド サンプルがない場合は、作物の集計統計情報を活用できます。世界の多くの地域では、作物の統計情報が集計され、定期的に公開されています。米国では、米国農務省の National Agricultural Statistics Service(NASS)が、郡ごと、主要作物ごとの詳細な作物統計情報を公開しています。2022 年のアイオワ州セロ ゴード郡のトウモロコシの作付面積は 161,500 エーカー、大豆の作付面積は 110,500 エーカーでした。

この情報から、上位 2 つのクラスタのうち、面積が最も大きいクラスタはトウモロコシ、もう 1 つはダイズである可能性が高いことがわかりました。これらのラベルを割り当てて、計算された面積と公開されている統計情報を比較してみましょう。

var cornFeature = selectedFc.sort('area', false).first();
var soybeanFeature = selectedFc.sort('area').first();
var cornCluster = cornFeature.get('cluster');
var soybeanCluster = soybeanFeature.get('cluster');

print('Corn Area (Detected)', cornFeature.getNumber('area').round());
print('Corn Area (From Crop Statistics)', 163500);

print('Soybean Area (Detected)', soybeanFeature.getNumber('area').round());
print('Soybean Area (From Crop Statistics)', 110500);

作物地図を作成する

これで、各クラスタのラベルがわかったので、各作物の種類についてピクセルを抽出し、それらを統合して最終的な作物地図を作成できます。

// Select the clusters to create the crop map
var corn = clustered.eq(ee.Number.parse(cornCluster));
var soybean = clustered.eq(ee.Number.parse(soybeanCluster));

var merged = corn.add(soybean.multiply(2));
var cropVis = {min: 0, max: 2, palette: ['#bdbdbd', '#ffd400', '#267300']};
Map.addLayer(merged.clip(geometry), cropVis, 'Crop Map (Detected)');

結果を解釈しやすくするために、UI 要素を使用して地図に凡例を作成して追加することもできます。

// Add a Legend
var legend = ui.Panel({
  layout: ui.Panel.Layout.Flow('horizontal'),
  style: {position: 'bottom-center', padding: '8px 15px'}});

var addItem = function(color, name) {
  var colorBox = ui.Label({
    style: {color: '#ffffff',
      backgroundColor: color,
      padding: '10px',
      margin: '0 4px 4px 0',
    }
  });
  var description = ui.Label({
    value: name,
    style: {
      margin: '0px 10px 0px 2px',
    }
  });
  return ui.Panel({
    widgets: [colorBox, description],
    layout: ui.Panel.Layout.Flow('horizontal')}
)};

var title = ui.Label({
  value: 'Legend',
  style: {fontWeight: 'bold',
    fontSize: '16px',
    margin: '0px 10px 0px 4px'}});

legend.add(title);
legend.add(addItem('#ffd400', 'Corn'));
legend.add(addItem('#267300', 'Soybean'));
legend.add(addItem('#bdbdbd', 'Other Crops'));


図: トウモロコシと大豆の作物が検出された作物地図

結果を確認する

衛星エンベディング データセットを使用し、地域の集計統計とローカル情報のみを使用して、フィールド ラベルなしで作物タイプ地図を取得できました。結果を USDA NASS Cropland データレイヤ(CDL)の公式の作物タイプ地図と比較してみましょう。

var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.date(startDate, endDate))
  .first();
var cropLandcover = cdl.select('cropland');
var cropMap = cropLandcover.updateMask(croplandMask).rename('crops');

// Original data has unique values for each crop ranging from 0 to 254
var cropClasses = ee.List.sequence(0, 254);
// We remap all values as following
// Crop     | Source Value | Target Value
// Corn     | 1            | 1
// Soybean  | 5            | 2
// All other| 0-255        | 0
var targetClasses = ee.List.repeat(0, 255).set(1, 1).set(5, 2);
var cropMapReclass = cropMap.remap(cropClasses, targetClasses).rename('crops');

var cropVis = {min: 0, max: 2, palette: ['#bdbdbd', '#ffd400', '#267300']};
Map.addLayer(cropMapReclass.clip(geometry), cropVis, 'Crop Landcover (CDL)');

Google の結果と公式地図には差異がありますが、最小限の労力でかなり良い結果が得られたことがわかります。結果に後処理ステップを適用することで、ノイズの一部を除去し、出力のギャップを埋めることができます。

Earth Engine コードエディタでこのチュートリアルの完全なスクリプトを試す