【実装編その１】50年ぶり！！スペインラ・パルマ島の溶岩流を衛星画像から捉えてみた

2022年2月3日 20:39

先日、スペイン領ラ・パルマ島で発生した火山噴火について、人工衛星の画像から溶岩流出範囲の抽出を行った。

溶岩流出範囲の抽出にはGoogle Earth Engineと呼ばれる人工衛星の統合プラットフォームを使って解析をしている。

Google Earth Engineは無料で人工衛星を利用することが可能で、プログラミングすることで人工衛星画像を加工・解析することも可能なプラットフォームだ。

本記事では、Google Earth Engine上で人工衛星のデータを取得し、解析結果を出力するまでの流れをまとめた。本記事の目標はGoogle Earth Engineを使用してスペイン領ラ・パルマ島の溶岩流出範囲の抽出ができるまでとする。

もし、人工衛星のデータを実際触れてみたい人は、ぜひこの記事を参考にして試してみてほしい。

1. Google Earth Engineの使い方

まずはGoogle Earth Engineにユーザー登録をして、Google Earth Engineの基本的な動作を習得しよう。

Google Earth Engineの基本的な使い方は以下の記事が分かりやすく載っているため、ぜひ参考してユーザー登録をしてみよう。

Google Earth EngineはGmailのメールアドレスとパスワードがあれば登録可能だ。記事に沿って、人工衛星画像を表示できればGoogle Earth Engineを使う準備が整う。

2. 溶岩検出の実装

続いて、溶岩流出範囲の抽出コードを示す。Google Earth EngineはJavaScriptでプログラミングすることで衛星データの加工・解析が可能だ。

今回は以下の記事を参考にして、溶岩流出範囲の抽出コードを作成した。

この記事では去年11月に沖縄に漂着した、軽石の抽出を試みている。
こちらも非常に興味深く参考になった記事なので、是非見ていただきたい。

溶岩流出範囲の抽出のソースは以下のリンクに格納されている。

上のリンクを開くと、図1のような画面が表示されるが、試しに何も考えずこのコードを動かしてみよう。画面中央タブの上部にある「Run」ボタンを押すと、プログラムが動き計算を開始する。
(なお、この画面を開くためには事前にユーザ登録が完了している必要がある)

計算が終わると、図2のような抽出結果が画面下部の地図に表示される。

抽出結果は以下の4つのLayerが表示される。
・clasified_after：噴火後の溶岩検出結果
・clasified2_after：噴火後のSentinel-2衛星画像
・clasified_before：噴火前の溶岩検出結果
・clasified2_before：噴火前のSentinel-2衛星画像

「Layers」タブにある各層のチェックボックスによって各Layerの表示を確認することが可能だ。

Layerを切り替えることで、噴火前後の比較や抽出結果の正誤を確認することが可能だ。

3. 検出コード

最後に、作成した溶岩検出コードを示す。
以下のコードが溶岩検出に使用したJavaScriptのコードだ。

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

var bands = ["B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8", "B8A", "B9", "WVP"];
// get mask image of landcover
var landMask = ee.Image("CGIAR/SRTM90_V4");

// get images defore eruption
var datasets_before = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate('2021-08-25', '2021-09-30')
                  .filterBounds(roi)
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
                  .map(maskS2clouds)
                  .map(function(image) {return image.updateMask(landMask)});
// get images after eruption
var datasets_after = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate('2022-01-01', '2022-01-15')
                  .filterBounds(roi)
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
                  .map(maskS2clouds)
                  .map(function(image) {return image.updateMask(landMask)});

// get first image
var dataset_before = ee.Image(datasets_before.first());
var dataset_after = ee.Image(datasets_after.first());

// classifiy
var traningPoints = cool_lava.merge(hot_lava).merge(mountain).merge(field).merge(city);

var traning = dataset_after.select(bands).sampleRegions({
 collection: traningPoints,
 properties: ['value'],
 scale: 10
})

var classifier = ee.Classifier.smileRandomForest(5).train({
 features: traning,
 classProperty: 'value',
 inputProperties: bands,
})

var classified_before = dataset_before.select(bands).classify(classifier);
var classified_after = dataset_after.select(bands).classify(classifier);

// set visualization
var visualization = {
  min: 0.0,
  max: 0.3,
  bands: ['B12', 'B11', 'B4'],
};

var classifiedVizParams = {
 min: 1,
 max: 6,
 // cool_lava, hot_lava, mountain, field, city
 palette: ['#ff6734', '#ff0000', '#00ff00', '#9aff55', '#f6f600'],
};

Map.addLayer(dataset_before.clip(roi), visualization, 'sentinel2_before');
Map.addLayer(classified_before.clip(roi), classifiedVizParams, "clasified_after");
Map.addLayer(dataset_after.clip(roi), visualization, 'sentinel2_after');
Map.addLayer(classified_after.clip(roi), classifiedVizParams, "clasified_after");

溶岩検出コードでは、噴火前後のSentinel-2衛星画像を取得して、機械学習で溶岩の検出を行った。機械学習にはランダムフォレストを用いている。

ただし、このコードをそのまま動かせば検出が可能になるわけではない。
検出の対象範囲の設定や機械学習のための教師データの選択など、コードを動かすためにやることがいくつかある。

そこで次回の記事ではこのコードの動かし方についてまとめていきたい。

この記事が気に入ったらサポートをしてみませんか？

【実装編 その１】50年ぶり！！スペイン ラ・パルマ島の溶岩流を衛星画像から捉えてみた

1. Google Earth Engineの使い方

2. 溶岩検出の実装

3. 検出コード

【実装編その１】50年ぶり！！スペインラ・パルマ島の溶岩流を衛星画像から捉えてみた