文系のための「正規化指標」（２）

リモートセンシングにおける解析手法として、

正規化指標は基礎的な手法の一つであり、

高度な分析においても部分的に使用されることもある。

NDVI（Normalized Difference Vegetation Index）がよく知られているが、

植生以外の状態を分析するための正規化指標もある。

そこで、今回は代表的な正規化指標を紹介することにする。

とりあえずは、使用するデータを準備する。

まず、検索バーからLandsat 8の Tier 1 Raw Scenenes を検索し、
インポートして変数名を landsat8としておく。

また、表示範囲を絞り込むために、
ジオメトリ・ツール使って、適当な場所にポイントを作成し、
作成した変数名を geomtry としておく。

var image = ee.Image(landsat8
    .filterBounds(geometry)
    .filterDate('2015-06-01', '2015-09-01')
    .sort('CLOUD_COVER')
    .first()
);

ついでに、表示用のカラーパレットも準備しておく。

var waterPalette = ['red', 'yellow', 'green', 'blue'];
var builtPalette = ['blue', 'green', 'yellow', 'red'];
var snowPalette = ['red', 'green', 'blue', 'white'];

まずは、植生の含水量の指標としてGao (1996) によって開発された

NDWI（Normalized Difference Water Index）という指標である。

この指標は近赤外線（NIR)と短波長赤外線（SWIR）の放射輝度を使用する。

NDWI = (NIR - SWIR) / (NIR + SWIR)

短波長赤外域は植物の葉の水分によって吸収されるという特性があるが、

近赤外域は植物の細胞を傷つけてしまうために反射される。

この特性によって植物の含水量を分析する方法がNDWIである。

一応、実際に適用する方法は以下の通り…。
ただし、NDWIはあくまでTMセンサを前提とした指標であり、
Landsat 8 に搭載されているOLIセンサの場合には要注意！

var ndwi = image.normalizedDifference(['B5', 'B6']);
Map.addLayer(
    ndwi, 
    {
        min: 0.0, 
        max: 0.3, 
        palette: waterPalette
    }, 
    'NDWI'
);

このコードの実行結果は以下の通り。



NDWBI（Normalized Difference Water Body Index）という指標もある。

この指標はNDWIと同じ1996年にMacFeeters (1996)によって発表された指標であり、

緑色光域と近赤外域の反射・吸収特性を使用して水域を抽出する。

NDWBI = (Green - NIR) / (Green + NIR)

水域は緑色光域においても緑色の光を反射するが、近赤外域は吸収する。

したがって、植生のある部分では近赤外域の反射が相対的に強くなり、

逆に、水域では近赤外の反射が相対的に弱くなる。

ということで、Google Earth Engine で実装してみる。

var ndwbi = image.normalizedDifference(['B3', 'B5']);
Map.addLayer(
    ndwbi, 
    {
        min: -0.5, 
        max: 0.2, 
        palette: waterPalette
    }, 
    'NDWBI'
);

そして、実行結果は以下の通り。



NDWIとNDWBIは名前もよく似ていて、両方とも水に関わる指標であるが、

前者は植物の葉に含まれる水分量、後者は水域を抽出するための指標である。

実際に両者を比較すると、どのような違いが見られるだろうか？

正規化指標は植物や水に関わる分析以外にも利用される。

そうした指標の一つがNDBI（Normalized Difference Built-up Index）である。

この指標は都市域を抽出するために用いられる指標である。

NDBI = (SWIR - NIR) / (SWIR + NIR)

この指標のモデルを見てみるとNDWIの式の関係が逆転していることが解る。

この指標は都市部の短波長赤外域の反射強度が近赤外域よりも高いことを利用している。

ということで…。

var ndbi = image.normalizedDifference(['B6', 'B5']);
Map.addLayer(
    ndbi, 
    {
        min: -0.4, 
        max: 0.1, 
        palette: builtPalette
    }, 
    'NDBI'
);

これも実行すると…。



地球温暖化に焦点を当てるのであれば、雪の分析が必要なるかもしれない。

その指標に利用されるのがNDSI（Normalized Difference Snow Index）である。

この指標は雪が短波長赤外域を吸収すること利用している。

NDSI = (Green - SWIR) / (Green + SWIR)

波長域とバンド名が解っていれば簡単であるが…。

var ndsi = image.normalizedDifference(['B3', 'B6']);
Map.addLayer(
    ndsi, 
    {
        min: -0.3, 
        max: 0.3, 
        palette: snowPalette
    }, 
    'NDSI'
);

そして、最後に…。

https://code.earthengine.google.com/5acf990cf5f81ad5c265c36bb07c5cf0

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

文系のための「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.

文系のための「GEEでのイメージの切り出し」

Google Earth Engine（GEE）における最大の魅力はその処理能力である。
全世界規模での解析を一瞬にして実現できてしまう。
これは従来のリモートセンシングとは大きく異なる点である。

とは言うものの…解析範囲は最低限で…ということもある。

GEEの処理速度は驚異的ではあるけれど、
無駄に広い範囲のデータを準備してしまうと、
さすがに、表示速度は遅くなってしまう。

ということで、今回は必要な範囲を切り出す方法を整理する。
実際には色々と方法はあるけれど、Google Driveへの保存方法と同様に、
ジオメトリの範囲で切り出す方法を実験してみることにする。

ジオメトリはジオメトリツールも使えるけれど、今回は手入力で準備する。

var area = ee.Geometry.Polygon(
    [
      [135.12,34.76],
      [135.12,34.6],
      [135.4,34.6],
      [135.4,34.76],
      [135.4,34.76],
      [135.12,34.76]
    ]
);

次に画像を準備する。今回も神戸市周辺のLandsat 8の画像にしておく。

今回は、Landsat 8 (Collection 1 Tier 1 TOA Refectance) のデータをインポートし、

変数名を「landsat8」としておく。詳細は以前の記事を参照…。

var ln8_2016 = ee.Image(

                landsat8

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

                .filterBounds(area)

                .sort("CLOUD_COVER")

                .first()

);

では、フィルタリングをかけて抽出した衛星画像をジオメトリで切り出す。

var ln8_2016_clip = ln8_2016.clip(area);

処理方法は非常に簡単である。Image オブジェクトの clip() オペレーションを使用し、

パラメータにジオメトリを指定するだけである。

今回は、ポリゴンとして作成した幾何オブジェクト「area」を指定している。

最後にこの結果を表示してみる。

Map.centerObject(area);

Map.addLayer(area);

Map.addLayer(ln8_2016_clip,

  {

    bands:["B4", "B3", "B2"],

    min: 0.0,

    max: 0.2

  });

以下が表示結果。無事に作成したオブジェクトによって切り出されている。


https://code.earthengine.google.com/d5813f0c376bf5d6f468833d075db1f1

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

文系のための「SENTINEL-2」

GEEでは様々なデータを検索し、使用することができる。
使いやすさの点では、やはり、Landsat シリーズであるが、
解像度が低いという問題がある。

そこで、今回はもう少し解像度が高い地球観測衛星を使ってみる。

実は、商用衛星を含めると解像度の高い人工衛星は色々と存在するが、
Google Earth Engine の中で使えるものに限ると、
Sentinel-2の衛星画像が最も使い勝手が良いように思える。

Sentinel-2は欧州宇宙機関（ESA: European Space Agency）という機関による
「コペルニクス計画」という地球観測研究を支えている人工衛星であり、
植生、地質、水域、海岸線などのモニタリングに用いられている。

Sentinel-2は2015年6月23日に打ち上げられた「Sentinel-2A」と、
2017年3月7日に打ち上げられた「Sentinel-2B」との二機で構成され、
5日周期でマルチスペクトル画像を取得することができる。

Sentinel-2によって得られたデータは13バンドのマルチスペクトルであり、
空間解像度はバンドごとに異なっているが、可視領域と近赤外については、
10m 解像度の画像を得ることができる。かなり使えそうである。

…ということで、早速、使ってみることにする。

GEEでこのセンサーで取得されたデータを使用するためには、
Code Editor の検索バーで「Sentinel-2」と検索し、RASTERS の項目から
現れた候補の中から「Sentinel-2 MSI: MultiSpectral Instrument, Level-1C」を選ぶ。



現れた候補をクリックすると、このデータセットの説明を見ることができる。
内容を確認できたら「Import」ボタンを押すと、
Code Editor に「imageCollection」という名前でインポートされる。

インポート名はそのままでも良いが、変数は解りやすい方が良い。
ということで、今回は規定の「imageCollection」という変数を
「sentinel_2」という変数名に変更しておくことにする。

では、さっそく、Sentinel-2の画像を扱ってみる。

Landsat の場合と同様に、最初に使用する画像を絞り込む。

今回はメタデータによるフィルタリングを試してみる。

var s2_filtered = ee.ImageCollection(sentinel_2

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

);

メタデータの内容を使ってフィルタリングを行うには、

filterMetadata(要素名, 演算子, 値) オペレーションを使用し、

条件に当てはまる画像のみを抽出することができる。

今回の場合、検索メタデータの要素名が「CLOUDY_PIXEL_PERCENTAGE」 で、

演算子が 「not_greater_than」、値が「5」となっているので、

「雲量比率が5%未満の画像のみを抽出せよ」という意味になる。

検索メタデータは扱う衛星画像（あるいは画像）によって異なるため、

データをインポートする前に説明を読み、

事前に使えそうなメタデータの要素を確認しておくことが重要である。

ちなみに、使用できる演算子は以下の通り…。

• "equals"（条件に等しい）
• "less_than"（条件より小さい）
• "greater_than"（条件より大きい）
• "not_equals"（条件に等しくない）
• "not_less_than"（条件より小さくない）
• "not_greater_than"（条件より大きくない）
• "starts_with"（条件を開始とする）
• "ends_with"（条件を終了とする）
• "not_starts_with"（条件を開始としない）
• "not_ends_with"（条件を終了としない）
• "contains"（条件を含む）
• "not_contains"（条件を含まない）

この辺りの処理の意味が解れば、メタデータの重要性が解るであろう。

上手く組み合わせることで、効率的に必要な画像を探し出し、

抽出された画像のみを使った分析を行うことができる。

今回の場合はシンプル・コンポジション画像を構築することにする。

Landsat の場合には特別の関数が用意されていたが、

Sentinel-2の画像を用いる際にはちょっとした工夫が必要となる。

とりあえず、以下のようにする。

var s2_filtered_med = s2_filtered.reduce(ee.Reducer.median());

ここでは、メタデータで抽出された ImageCollection （s2_filtered）に対し、

reduce()というオペレーションを実行している。

実行している内容は ee.Reducer.median() である。

要するに、ImageCollection として集められた画像から中央値（Median）を抽出する。

もちろん、中央値以外にも様々な演算子が存在しているが、

参考までによく使いそうな演算子を挙げておく。

• ee.Reducer.sum()：合計値を返す。
• ee.Reducer.mean()：算術平均を返す。
• ee.Reducer.median()：中央値を返す。
• ee.Reducer.min()：最小値を返す。
• ee.Reducer.max()：最大値を返す。
• ee.Reducer.minMax()：最小値と最大値を返す。
• ee.Reducer.variance()：分散を返す。
• ee.Reducer.stdDev()：標準偏差を返す。
• ee.Reducer.count()：Nullでない値の値の数を返す。
• ee.Reducer.first()：最初の入力の値を返す。
• ee.Reducer.last()：最後の入力の値を返す。

さて、この処理の結果をConsoleで確認してみる。

print(s2_filtered_med);

結果をみてみると、各バンドの名前の後ろに「_median」とついている。

Reducerを使って処理を行うと、バンド名の後ろに演算子の名前がつく。

Reducerを実行したい場合には Console に出力して確認することが重要である。

そして、最後にSentinel-2をマップ上に表示する。

Map.addLayer(s2_filtered_med,

  {

    bands: ['B4_median', 'B3_median', 'B2_median'], 

    min: 0,

    max: 2500

  },'Composite_Med'

);

最終的に表示された画像をみてみると、確かにLandsatの画像よりも細かい。

https://code.earthengine.google.com/f0e45f6d71181e0e2a8f1ef4664c7064

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

文系のための「ALOS/AVNIR-2」

地球観測衛星の保有国はアメリカやヨーロッパだけではない。
日本も複数の地球観測衛星を打ち上げていて、
国内外の研究で活用されている。

ということで…

今回は国産の人工衛星である「だいち（ALOS）」のデータを使ってみる。
この人工衛星にはPRISM、AVNIR-2、PALSARの三つの観測機器が搭載されていて、
今回はAVNIR-2（高性能可視近赤外放射計2型）に焦点を当ててみる。

AVNIR-2センサは1996年に打ち上げられた「みどり」に搭載されたセンサー、
AVNIRの改良型であり、RGBの可視光に加えて近赤外の波長帯を観測する。
AVNIR-2の基本スペックは以下の通り…。

• 観測波長帯：
• Band 1 : 0.42 〜 0.50μm（青色光）
• Band 2 : 0.52 〜 0.60μm（緑色光）
• Band 3 : 0.61 〜 0.69μm（赤色光）
• Band 4 : 0.76 〜 0.89μm（近赤外光）
• 地上分解能：
• 10m
• 観測幅：
• 70km

Landsat 8 と比較すると波長帯の数が少ないが、

地上分解能が高いため、解像度の高い画像を得ることができ、

土地被覆分類図や土地利用分類図などへの利用が想定されている。

GEEでこのセンサーで取得されたデータを使用するためには、

Code Editor の検索バーで「ALOS/AVNIR-2 ORI」と検索する。

すると、RASTERS の項目に候補が現れる。

現れた候補をクリックすると、このデータセットの説明を見ることができる。

説明をみてみると、各バンドの値は1〜255に正規化されていることが解る。

内容を確認できたら「Import」ボタンを押して、インポートする。

インポート名はそのままでも良いが、変数は解りやすい方が良い。

ということで、今回は規定の「imageCollection」という変数を

「avnir2」という変数名に変更しておくことにする。

Landsat の場合には「Cloud Score」のメタ情報を上手く使って、

雲量が少ない順にデータを並べ替えることができたが、

ALOSにはそうした情報がメタデータに含まれていない…。

そのため、必要なデータに関しては何らかの方法で検索する必要がある。

今回はJAXAのALOSのデータセット検索ツールを使用する。

とりあえず、以下のページに移動してみる…。

http://www.eorc.jaxa.jp/ALOS/alos-ori/index.html

「ALOS オルソ補正画像プロダクト(ALOS-ORI)」と書かれた直下に

Google Maps が表示されていて、その地図上に観測範囲と

その中心にピンが配置されているハズ。

そして、必要な範囲に相当する場所のピンをクリックすると、

その範囲を撮影した範囲のリストがポップアップで表示されるので、

最上段の「View all」をクリックしてリストを開く。

リストを開くと、サムネイル付きで衛星画像の詳細を確認できるので、

このリストを確認しながら雲量の少ない画像を探す。

そして、使えそうな画像を発見したら「SCENE ID」を控えておく。

とりあえず、近畿地方の三つの画像をピックアップしてみた…。

• ALAV2A172272900
• ALAV2A122822900
• ALAV2A127052900

本来であれば近い時期のものを選ぶ必要があるが、

雲が出やすい地域では、少々、厄介である。

ということで、今回は練習なので気にしないことにする。

再び、Code Editor に戻って、次のように変数を作成する。

var kinki_1 = ee.Image(avnir2
  .filterMetadata('SCENE_ID', 'equals', 'ALAV2A172272900')
  .first()
);

var kinki_2 = ee.Image(avnir2
  .filterMetadata('SCENE_ID', 'equals', 'ALAV2A122822900')
  .first()
);

var kinki_3 = ee.Image(avnir2
  .filterMetadata('SCENE_ID', 'equals', 'ALAV2A127052900')
  .first()
);

ここでポイントとなるのは .filterMetadata 関数の部分である。
メタデータの「SCENE_ID」が一致（equals）となるように、
フィルタリングを行って必要なImageCollectionを抽出する。

ただし、この処理によって検索されるハズの画像は一枚だけなので、
最後に .first() 関数を使ってImageCollection からImageに変換している。
少し、ややこしいかもしれない。

これで必要な画像を表示することができるが、分析するには不便である。
そのため、三つの Image を新しい ImageCollection として統合し、
そこからモザイク画像を作成してみる。

// Create a ImageCollection by using extracted three images.
var kinki = ee.ImageCollection([kinki_1,kinki_2,kinki_3]);

// Create a mosaic from the ImageCollection.
var kinki_mosaic = kinki.mosaic();

// Add the result to the map layer.
Map.addLayer(
  kinki_mosaic,
    {
      bands: ['B4', 'B3', 'B2'],
      max: 255,
      min: 1
    }
);

これを実行すると以下のようになる。

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

文系のための「GEEでのイメージの保存」

Google Earth Engine（GEE）で作成した画像や、
分析結果を他のソフトウェアで利用するにはどうすれば良いのか？
おそらく、多くの人が感心を持つ問題であろう。

GEEで画像を出力する方法は他の一般的なソフトウェアとは大きく異なる。
というのは、作成した画像を直接的にダウンロードすることはできず、
一旦、クラウド・ストレージ上に出力してから、ダウンロードしなけれなならない。

GEEからの出力先としては三種類が用意されている。

• Google Drive
• Google Storage Cloud
• Google Earth Engine Assets

一般ユーザーの多くはGoogle Driveへの出力を行うことになるだろう。

まずは、検索バーを使って、Landsat 8 のImageCollectionをインポートし、
インポートしたImageCollectionの名前を「landsat8」に変更しておく。
そして、時期と場所でフィルタリングし、雲がもっとも少ない一枚を取得する。

その次に画像の抽出範囲および保存範囲を既定するためのジオメトリを作成する。
とりあえず、神戸市周辺のエリアを抽出してみる。
ジオメトリ・ツールを使ってインポートすることもできるが、今回は手打ち入力。

var area = ee.Geometry.Polygon(
    [
      [135.12,34.76],
      [135.12,34.6],
      [135.4,34.6],
      [135.4,34.76],
      [135.4,34.76],
      [135.12,34.76]
    ]
);

この部分に関しては文系のための「GEEの幾何オブジェクト（Geometry）」を参照。

次に、インポートしたImageCollectionから必要なデータを抽出する。

var ln8_2016 = ee.Image(
                landsat8
                .filterDate("2016-01-01","2016-12-31")
                .filterBounds(area)
                .sort("CLOUD_COVER")
                .first()
                );
Map.centerObject(ln8_2016, 10);

詳しい手順は文系のための「Image CollectionからのImageの抽出」を参照。

では、さっそく、抽出した画像をGoogle Driveへと出力する。
そのままでもマルチバンド・画像として出力することができるが、
今回はトゥルーカラー画像の状態で出力してみる。

var output_image = ln8_2016.select(['B4', 'B3', 'B2']);

最後に出力を行う。出力する方法は以下の通り。
出力するためには、生成したイメージ、範囲、出力解像度の指定が必要となる。
今回は１ピクセル30メートル解像度を指定する。

Export.image.toDrive({
            image: output_image,
            description: 'output',
            scale: 30,
            region: area
});

これまでに扱ったGEEの処理と異なる点が一つある。
出力は「タスク（Task）」として処理されるため、
タスク・マネージャー上から実行しなければならない。

タスク・マネージャーはCode Editor上段の右側のタブにある。



この画面の最上段にある「output」というのが登録されたタスクである。
右横にある「RUN」というボタンをクリックするとタスクが実行される。
実行されたタスクはクラウド上で処理されるため、ログアウトしても問題ない。

タスクを実行すると以下のようなダイアログが出るので必要があれば修正する。



そして、「Run」をクリックするとタスクが実行される。
範囲が広いと長い時間がかかることがあるが、根気強く待つ…。
無事に終了すると、Google Driveの指定されたフォルダに出力される。

Google Drive に保存されたファイルは位置情報を持ったGeoTIFFという画像であるので、
たとえば、QGISなどのGISソフトで開くと正しい位置に読み込まれる。
これで必要に応じてその他の分析とも組み合わすことができるだろう。