文系のための「GEEの外部モジュールの作成と利用」

Google Earth Engine（GEE）は慣れると使いやすいし、心地よい。
若い頃は新しいソフトウェアの操作をすぐに覚えたけれど、
歳をとると、操作方法を覚えるのは一苦労である…。

その点、RやGEEは操作方法を覚える必要が無いので良い。

しかしながら、慣れれば、慣れるほど、
同じコードを書くのが面倒になってくる。
できれば、コードの再利用をしてみたくなる。

例えば、カラーパレットの定義やUIなどを毎回作るのは面倒である。

そこで、今回はよく使う機能をモジュール化し、
再利用するための方法を検討する。
上手く使いこなすことができれば生産性も向上すると期待できる。

まず、分析用のモジュールを格納するためのレポジトリを用意する。
Scripts タブの NEW をクリックし、オブジェクトの新規 メニューを開き、
Repository を選択する。



レポジトリの名前は何でも良いけれど、今回は「modules」にしておく。



すると、Scripts タブに新しいレポジトリが登録される。
今回はこのレポジトリの中に二つのプログラムを登録することにする。
まずは、可視化パラメータを作ってみる。

新規ファイルを作成し、以下の内容を貼り付ける。
コードはカラーパレットやLandsatの可視化パラメータである。
貼り付けたら「visParam」という名前で作ったレポジトリに保存する。

exports.veg = ['blue', 'red', 'yellow', 'green'];
exports.dem = ['green', 'yellow', 'orange', 'brown', 'white'];

exports.rb = ['red', 'white', 'blue'];
exports.br = ['blue', 'white', 'red'];
exports.ro = ['red', 'yellow', 'orange'];

exports.Ln5TOA = {
  bands: ['B3_median', 'B2_median', 'B1_median'],
  gamma: 1,
  max: 0.2,
  min: 0.0
};

exports.Ln8TOA = {
  bands: ['B4_median', 'B3_median', 'B2_median'],
  gamma: 1,
  max: 0.2,
  min: 0.0
};

これまでのは「var」によって変数を宣言していたが、
上記のコードでは変数の前に「exports.」が付いている。
この魔法の言葉をつけることで外部から関数や変数を呼び出すことできる。

ついでに、MapLinkerを作成するコードについてもモジュール化してみる。
新規ファイルを作成、以下のコードを貼り付ける。
中身の処理については以前の内容を参照。

今回作成する関数「mapLinker()」は左右に画像を並べて表示するプログラムであり、
画像オブジェクト（left と right）、可視化パラメータ（l_visParam と r_visParam）、
および、表示用のラベル名（l_label と r_label）の6つを引数とする。

exports.mapLinker = function(
　　　　left, right, l_visParam, r_visParam, l_label, r_label) {
    // 複数のマップオブジェクトを格納するための空の配列を定義する。
    var maps = [];
   
    // マップ領域を準備する。
    var map_left = ui.Map();        // 左側のマップ領域
    var map_right = ui.Map();       // 右側のマップ領域
   
    // それぞれのレイヤにそれぞれの画像を追加する。
    map_left.addLayer(left, l_visParam, l_label);
    map_right.addLayer(right, r_visParam, r_label);   

    // それぞれのMapにタイトルラベルを追加する。
    map_left.add(ui.Label(l_label));
    map_right.add(ui.Label(r_label));
   
    // それぞれのMapにレイヤー・マネージャーを追加する。
    map_left.setControlVisibility(true);
    map_right.setControlVisibility(true);   
   
    // それぞれのMapを配列に追加する。
    maps.push(map_left);
    maps.push(map_right);

    // マップ・リンカーオブジェクトを作って、複数のマップを連動させる
    var linker = ui.Map.Linker(maps);
   
    // UIをリセットする。
    ui.root.widgets().reset(maps);

    // 地図の中心をひとつ目のマップのオブジェクトの衛星画像に設定する。
    maps[0].centerObject(map_left, 9);
};

このように「customUi」と「visParam」の二つのモジュールを準備する。

これらのモジュールの準備ができたら、次に、これらを実行するプログラムを準備する。

とりあえず、Landsat 5 と Landsat 8 の画像（Tier 1 TOA Reflectance）を検索し、

それぞれ、「landsat5」と「landsat8」という名前でインポートする。

その上で、以下のコードを貼りけて実行する。

// 外部ライブラリを読み込む
var customUi = require('users/yufujimoto/modules:customUi');
var visParam = require('users/yufujimoto/modules:visParam');

// Landsat 5 の画像を読み込む。

var ln5 = ee.ImageCollection(landsat5

  .filterDate("1989-01-01","1990-12-31")

  .filterMetadata("CLOUD_COVER","not_greater_than",10)

);

// Landsat 8 の画像を読み込む。

var ln8 = ee.ImageCollection(landsat8

  .filterDate("2016-01-01","2017-12-31")

  .filterMetadata("CLOUD_COVER","not_greater_than",10)

);

// シンプルコンポジットを作成する。

var ln5_sc = ln5.reduce(ee.Reducer.median());

var ln8_sc = ln8.reduce(ee.Reducer.median());

// 作成した外部モジュールを使ってMapLinkerを作成する。
customUi.mapLinker(
        ln5_sc, ln8_sc, 
        visParam.Ln5TOA, visParam.Ln8TOA, 
        "1989-1990", "2016-2017"
);

ここで重要となるのが最初の二行と最後の一行の部分である。
すなわち、require() 関数によって読み込んでインスタンス化し、

インスタン化されたオブジェクトを通して実行する。

すなわち、可視化パラメータに関しては visParamという名前でインスタンス化し、
visParam.Ln5TOA と visParam.Ln8TOA としてモジュールで定義した変数を呼び出す。
同様に、customUi.mapLinker() によって独自に定義した関数を呼び出している。

実行した結果は以下の通り。モジュール化によってコードが見やすくなった。


https://code.earthengine.google.com/7af5f2f2c055afdb383d3ec9f297bf38



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

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

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.