文系のための「地すべり分析」（西日本豪雨災害）

2018年7月に発生した豪雨災害は西日本に甚大な被害を及ぼした。

この災害では多くの尊い命が失われただけでなく、

様々な形で被災地に大きな被害を与えることになった。

ちょうど、この時、私は鹿児島県三島村で調査をしていて、

調査終了後、新幹線の車窓から被害状況を目の当たりにし、

非常に大きなショックを受けたことを覚えている。

さて、こうした災害が発生すると、多くの人は最も大きな被災地に目を向ける。

もちろん、これは当然のことではあるけれど、そうした地域では手が回らず、

被害の実態を把握し、復旧するまでに非常に長い時間がかかることがある。

例えば、高知県大豊町を通る高知自動車道が土砂災害によって寸断され、

現在も復旧作業が続けられているが、未だに対面通行の状態が続き、

少なからず、陸運に影響を及ぼしている。

道路に限らず、大規模な災害が発生した際に被害状況を迅速に把握し、

復興計画および復興のための予算措置を考えることは極めて重要である。

しかしながら、この段階には非常に長い時間を要することになる。

…ということが、2018年12月26日にGoogle Japan で開催された Remo-thon において

「災害チーム」として集まったメンバーで考えた課題の一つとなった。

ちなみに、Remo-thon = Remote Sensing + Hacking + marathon という意味...。

話がそれてしまったが、とりあえず、我々のチームでは以下のことを考えた。

・地すべりの範囲を定量的に算出したい。

・解像度は高い方が望ましいため、Sentinel-2 を使用するのが妥当である。

・森林地域なのでNDVI（植生指標）を使える可能性が高い。

まず、「Sentinel-2 MSI: MultiSpectral Instrument, Level-1C」を検索する。
Sentinel-2の詳細については既に解説済みなので、そちらを参照。

とにかく、Code Editor にインポートした後、変数を「s2」として話を進める。

最初は、パラメータ類の設定を行う。これらは表示の時に使用する。

// NDVI 用の表示パラメータの準備
var visparam_ndvi = [
                     'FFFFFF', 'CE7E45', 'DF923D', 'F1B555',
                     'FCD163', '99B718', '74A901', '66A000',
                     '529400', '3E8601', '207401', '056201',
                     '004C00', '023B01', '012E01', '011D01',
                     '011301'
];

// 結果表示用の表示パラメータの準備

var visparam_diff = ['FFFFFF', 'FF0000'];

次に、分析範囲を定義する。ジオメトリ・ツールを使っても良いが、

今回は、手入力で範囲を設定する。こちらの方が都合が良いこともある。
分析範囲の切り出し方法についての詳細は既に説明している。

// 高知県大豊町の土砂災害地域の範囲

var area = ee.Geometry.Polygon(

    [

      [133.64689749333183, 33.83090489745764],

      [133.66440695378105, 33.83090489745764],

      [133.66440695378105, 33.84017316151725],

      [133.64689749333183, 33.84017316151725],

      [133.64689749333183, 33.83090489745764]

    ]

);

解析対象となるSentinel-2の画像を二種類準備する。

一方は、災害発生後の画像（2018年7月10日〜2018年12月25日）であり、

もう一方は、災害発生前の画像（2017年1月1日〜2018年12月31日）である。

// 災害後の画像

var s2_af = ee.ImageCollection(s2

  .filterBounds(area)

  .filterDate('2018-07-10', '2018-12-25')

  .filterMetadata("CLOUDY_PIXEL_PERCENTAGE","not_greater_than",5)

);

// 災害前の画像

var s2_nc = ee.ImageCollection(s2

  .filterBounds(area)

  .filterDate('2017-01-01', '2017-12-31')

  .filterMetadata("CLOUDY_PIXEL_PERCENTAGE","not_greater_than",5)

);

災害後の画像は「s2-af」という変数に格納され、

災害前の画像は「s2-nc」という変数に格納されている。

この段階でのデータ型は ee.ImageCollection となっている。

いずれの画像についても、雲の影響を最小限に抑えるために、

メタデータの情報を使ってフィルタリングを行い雲量5%以下を抽出しているが、

この段階では雲がかかった画像が含まれる可能性がある。

そこで、複数時期の画像のImageCollectionからシンプル・コンポジットを構築し、

さらに、最初に定義した分析範囲で画像を切り出す。

色々と方法はあるが、今回はシンプルに中央値（Median）を使用する。

// ImageCollectionにおける各ピクセルを中央値で集約する。

var s2_af_sc = s2_af.reduce(ee.Reducer.median()).clip(area);

var s2_nc_sc = s2_nc.reduce(ee.Reducer.median()).clip(area);

災害後の雲なし画像が「s2_af_sc」であり、

災害前の雲なし画像が「s2_nc_sc」となる。

データ型はReduceされてee.Imageになる。

// 災害後の画像 NDVI の算出を行う。

var s2_af_sc_ndvi = s2_af_sc.normalizedDifference(
                    ['B8_median', 'B4_median']
);

// 災害前の画像 NDVI の算出を行う。
var s2_nc_sc_ndvi = s2_nc_sc.normalizedDifference(
                    ['B8_median', 'B4_median']
);

災害後のNDVIと災害前のNDVIとの差分をとり、「diff_ndvi」に格納する。

// 災害後の NDVI から災害前の NDVI を引く。

var diff_ndvi = s2_af_sc_ndvi.subtract(s2_nc_sc_ndvi);

地すべりが発生した部分では地面が露出するため植生は無い状態になる。

したがって、災害後に地面が露出した部分は通常の状態よりも植生指標は低い。

すなわち、この計算によって得られた差分画像から値の低い場所を探せば良い。

問題は、災害の前後での変化量のしきい値の設定であるが、

今回は変化量が正規分布に従うと仮定した上で、

標準偏差を計算し、−2σ以下を地滑りエリアとすることにする。

標準偏差や算術平均などの計算を行うためにはReducerを利用し、

ピクセル値の集計については、ImageオブジェクトのReduceRegionを用いる。

今回は標準偏差と算術平均を使用し、それぞれの値を計算する。

ReduceRegion は複数のバンドに対して適用され、その返り値は配列となる。

したがって、得られた結果を取り出すためには、get() で取り出す必要がある。

NDVIの差分画像の場合はバンド名「nd」を使って、 get("nd") で取り出す。

なお、数値は ee.Number() を使って型変換を行う必要がある。

// 差分データの標準偏差を計算する。

var diff_ndvi_sd = ee.Number(

  diff_ndvi.reduceRegion({

    reducer: ee.Reducer.stdDev(),

    geometry: area,

    scale: 10,

  }).get("nd")

);

// 差分データの算術平均を計算する。

var diff_ndvi_mean = ee.Number(

  diff_ndvi.reduceRegion({

    reducer: ee.Reducer.mean(),

    geometry: area,

    scale: 10,

  }).get("nd")

);

ここまでの処理で必要な値は算出できたので、−2σでしきい値を設定する。

// 平均から標準偏差の2倍値を引く（−2σ）

var thresh = diff_ndvi_mean.subtract(diff_ndvi_sd.multiply(2));

// 計算されたしきい値を使ってNDVIの差分画像を2値化する

var result = diff_ndvi.lt(thresh);

最後に表示を行う。

// 分析範囲にズームする

Map.centerObject(area);

// レイヤに画像を追加する

Map.addLayer(
    s2_af_sc_ndvi, 
    {min: 0, max: 1, palette: visparam_ndvi}, 'NDVI'
);

Map.addLayer(
    s2_nc_sc_ndvi, 
    {min: 0, max: 1, palette: visparam_ndvi}, 'NDVI'
);

Map.addLayer(
    result, 
    {min: 0, max: 1, palette: visparam_diff}, 'Subtruct'
);

もしも、この結果をGISで利用するには、

以下のコマンドを追加し、タスクを実行すると、

Google Drive に出力結果を保存することができる。

// 結果をGoogle Driveに保存する

Export.image.toDrive({

  image: result,

  description: 'landslide area extraction',

  fileNamePrefix: 'result',

  region: area,

  scale: 10,

  crs: 'EPSG:3857'

});

https://code.earthengine.google.com/2c0a8803e722d27da8fe8e7c78fdc8c8

This work is licensed under a Creative Commons Attribution 4.0 International License.