叔母の旦那さんに、パキスタン出身の人がいる。

この前、パキスタン人の叔母の旦那さんと会話をしたら、
両親はパキスタンに住んでいるとのことだった。

今回のパキスタンの大洪水は、意外と身近な人が被害を受けてるなと感じるね。

先日、パキスタンが「歴史的な洪水」に見舞われた。
世界屈指の大河であるインダス川が大規模に氾濫し、「国土の3分の1が水没」したとも伝えられる。

この大規模な洪水は、人工衛星のデータからも、その規模の大きさが把握できた。

今回はパキスタンの大洪水について、人工衛星からの解析を試みた。
なお、今回の解析は、無料で人工衛星画像を解析できるGoogle Earth Engineを使用する。

1. 光学衛星を使用した洪水前後の比較

まず、人工衛星からパキスタンの洪水がどのように観察できるかをみてみよう。

以下は、パキスタンの中心を流れる、インダス川流域の人工衛星画像だ。
左は洪水発生前の画像であり、右は洪水発生後の画像だ。

使用した人工衛星は、Sentinel-2衛星と呼ばれる衛星だ。

Sentinel-2衛星は、欧州の宇宙機関であるESA(European Space Agency)が運用する光学衛星だ。Sentinel-2衛星のデータは無償で使用することが可能で、今回のような洪水検知にも活用することができる。

衛星画像を確認すると、インダス川流域の多くの箇所が洪水の被害を受けている様子が確認できる。

2.SAR衛星を使用した洪水前後の比較

次に、別の衛星でも、洪水の様子を確認してみよう。
ESAは「Sentinel-2衛星」以外にもいくつかの人工衛星を運用している。

その中で、洪水などの被害状況の把握に活用されるのが、「Sentinel-1衛星」だ。

「Sentinel-1衛星」は、合成開口レーダー(SAR)と呼ばれる特殊なレーダーを搭載した人工衛星だ。

SAR衛星の特徴は、天候などの影響を受けず、洪水の被害状況を観測できる点にある。

以下に「Sentinel-1衛星」で観測した洪水前後のパキスタンの画像を示す。
(撮影場所は「Sentinel-2衛星」の画像と同じ)

洪水発生前後のSAR画像を比較すると、洪水発生後に画像の「黒い箇所」が拡大していることが確認できる。

細かなSAR画像の説明は省略するが、SAR画像内の「黒い箇所」は洪水での浸水地域を示している。

洪水前後の比較により、撮影箇所一帯が洪水で水没したことがわかる。

3. パキスタンにおける洪水検出

2つの人工衛星の画像から、パキスタンにおける大規模な洪水の状況を確認することができた。

それでは最後に、この洪水の被害範囲の推定に、人工衛星がどのように活用されるのかを紹介しよう。

洪水において、人工衛星は「被害範囲の把握」に活用される。
人工衛星は、広範囲の情報を一度の撮影で把握できることが、強みの１つだ。

実際にそのメリットを活用して、広範囲に発生した洪水の被害状況を検出される。

以下の画像は、実際に衛星画像を使って、自動的に検出した洪水の被害地域を示す。

洪水被害地域の推定には、Sentinel-1衛星のSAR画像を用いて行った。
画像内のオレンジの箇所は洪水によって被害を受けた場所を示す。

人工衛星の画像を活用することで、洪水の被害範囲が自動的に検出することが可能だ。

<実行コード>
Google Earth Engine上で実行可能なコード

// Sentinel-2用の表示パラメータ準備
var visualization = {
min: 0.0,
max: 4000.0,
bands: ['B4_median', 'B3_median', 'B2_median'],
};

// 結果表示用の表示パラメータの準備
var visParams_diff = {min: 0,
                      max: 1,
                      bands: ['VV_median'],
                      palette:['FFFFFF', 'FF0000']};

// パキスタンの土砂災害地域の範囲
var area = ee.Geometry.Polygon(
    [
      [67.90287383594544,26.96842971392122],
      [67.93291457691224,26.96842971392122],
      [67.93291457691224,26.991070596923866],
      [67.90287383594544,26.991070596923866],
      [67.90287383594544,26.96842971392122]
    ]
);

// 災害前のSentinel-1画像
var s1_before = ee.ImageCollection("COPERNICUS/S1_GRD").filterBounds(area)
                      .filterDate('2022-08-02', '2022-08-06');
// 災害前のSentinel-2画像
var s2_before = ee.ImageCollection("COPERNICUS/S2_SR").filterBounds(area)
                      .filterDate('2022-08-03', '2022-08-05')
                      .filterMetadata("CLOUDY_PIXEL_PERCENTAGE","not_greater_than",50);

// 災害後のSentinel-1画像
var s1_after = ee.ImageCollection("COPERNICUS/S1_GRD").filterBounds(area)
                      .filterDate('2022-08-25', '2022-08-30');
// 災害後のSentinel-2画像
var s2_after = ee.ImageCollection("COPERNICUS/S2_SR").filterBounds(area)
                      .filterDate('2022-08-28', '2022-08-30')
                      .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'not_greater_than', 50);

// ImageCollectionにおける各ピクセルを中央値で集約
var s1_after_cent = s1_after.reduce(ee.Reducer.median()).clip(area);
var s1_before_cent = s1_before.reduce(ee.Reducer.median()).clip(area);
var s2_after_cent = s2_after.reduce(ee.Reducer.median()).clip(area);
var s2_before_cent = s2_before.reduce(ee.Reducer.median()).clip(area);

// 災害後のNDVIから災害前のNDVIを引く
var diff_s1 = s1_after_cent.subtract(s1_before_cent);

// 差分データの標準偏差を計算
var diff_sd = ee.Number(
  diff_s1.reduceRegion({
    reducer: ee.Reducer.stdDev(),
    geometry: area,
    scale: 10
  }).get('VV_median')
);

// 差分データの算術平均を計算する
var diff_mean = ee.Number(
  diff_s1.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: area,
    scale: 10
  }).get('VV_median')
);

// 平均から標準偏差を引く(-σ)
var thresh = diff_mean.subtract(diff_sd.multiply(1));

// 計算されたしきい値を使ってNDVIの差分画像を2値化する
var result = diff_s1.lt(thresh);

// 分析範囲にズームする
Map.centerObject(area)

// レイヤに画像を追加する
Map.addLayer(s2_after_cent, visualization, 'sentinel2_after');
Map.addLayer(s2_before_cent, visualization, 'sentinel2_before');
Map.addLayer(s1_before_cent, {min: -25, max: 20}, 'sentinel1_before');
Map.addLayer(s1_after_cent, {min: -25, max: 20}, 'sentinel1_after');
Map.addLayer(result, visParams_diff, 'Subtruct');