はじめに

Image編；
https://qiita.com/iwasaki_kenichi/items/0f8087b870688916525e

Google Earth Engineというものがあります。https://earthengine.google.com

Google Earth Engineは、数ペタバイト規模の衛星画像と地理空間データセットを惑星規模の解析機能と組み合わせて、科学者、研究者、開発者が変化を検出し、傾向をマッピングし、地球表面の違いを定量化できるようにします。 webサイトより引用

非常に興味深いので、使えるようになりたいと思っています。そこでGoogle Earth Engine（以降GEE）を使えるようになるため、Google が用意しているサンプルコードを写経していきたいと思います。ただし、ただの写経だとあんまり味気ないため、下記3つを自分に課したいと思います。
■Phase１：まずは普通に写経（コピペは厳禁）
■Phase２：コードの解説

※Image編では、
■Phase３：写経したコードを元に少し改良する（余裕があったら）
を入れていましたが、外します。ちょっと時間がかかる...

QiitaはGoogle Earth Engineのシンタックスハイライトに対応していないので、少し見にくいです。

サンプルコードについて

Google が用意しているサンプルコードは
- Image
- Image Collection
- Feature Collection
- Charts
- Arrays
- Primitive
- Cloud Masking
- Code Editor
- Demos
- User Interface
の10コースあります。

今回は、Image Collectionのサンプルコードを読み解いて行きたいと思います。

Image編は下記；
https://qiita.com/iwasaki_kenichi/items/0f8087b870688916525e

ImageCollection

ImageCollectionは、二番目のコースです。
ee.Imageは単一のデータを指定しましたが、ee.ImageCollectionはスタッキングされた大量のデータを扱います。

Clipped Composite

ImageCollectionにて取得した画像を使い、ee.FeatureCollectionにてフィルターし、
clip.ToCollectionにてフィルタしたデータをクリップします。
ここではFusion Tableというものを使っています。
Fusion TableとはGoogleが提供しているデータベースサービスのことです。
参考サイト：http://blog.thetheorier.com/entry/fusion-tables

Phase1

下記コードです。

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1')
        .filterDate('2000-04-01','2000-07-01')

var median = collection.median()

var fc = ee.FeatureCollection('ft:1fRY18cjsHzDgGiJiS2nnpUU3v9JPDc2HNaR7Xk8')
        .filter(ee.Filter.or(
          ee.Filter.eq('Name','Nevada'),
          ee.Filter.eq('Name','Arizona')));

var clipped = median.clipToCollection(fc)

Map.setCenter(-110,40,5)
var vis_Params= {band:['B3','B2','B1'],gain:[1.4,1.4,1.1]}

print(fc)

//Map.addLayer(fc,'fc')
Map.addLayer(clipped,vis_Params,'clipped composite')

Phase2

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1')
        .filterDate('2000-04-01','2000-07-01')

ee.ImageCollectionにて、Landsat7の2000年4月1日から2000年7月1日までのデータを読み込んでいます。
この状態での画像は以下の通りです。

var median = collection.median()

ここでは、Collection全体で各ピクセルのすべての値の中央値を計算し、Collectionを減らしています。
以下のような画像ができます。

中央値を採用しているので、極端に暗いところもなければ極端に明るいところもないような画像ができているかと思います。

var fc = ee.FeatureCollection('ft:1fRY18cjsHzDgGiJiS2nnpUU3v9JPDc2HNaR7Xk8')
        .filter(ee.Filter.or(
          ee.Filter.eq('Name','Nevada'),
          ee.Filter.eq('Name','Arizona')));

ee.FeatureCollectionにてGoogle Fusion Tableを呼び出しています。ee.FeatureCollectionはフィルタリングやソート、レンダリングなど、セット全体に対する操作を可能とします。
今回サンプルコードでどのようなFusion Tableなのかわかりませんが、後段処理で
ee.Filter.eq('Name','Nevada')やee.Filter.eq('Name','Arizona')など、州名でフィルタリングしているため、米国の全州都の名前があると考えられます。

Map.addLayer(clipped,vis_Params,'clipped composite')

最後に表示させます。

アリゾナとネバダの部分にクリップされた様子がわかります。

Expression Map

NDVIとSAVIを計算し、Map.addLayerにて表示させています。表示させる画像はee.ImageCollectionを使用しているため、スタックされた複数の画像を使用しています。
そのため、広範囲にわたったNDVIおよびSAVI画像が表示されています。

Phase1

以下コードです。

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_TOA')
      .filterDate('2002-11-01','2002-12-01')

var NDVI = function(image){
  return image.expression('float(b("B4") - b("B3")) / (b("B4") + b("B3"))')
}

var SAVI = function(image){
  return image.expression(
    '(1 + L) * float(nir - red) / (nir + red + L)',
    {
      'nir':image.select('B4'),
      'red':image.select('B3'),
      'L':0.2
    })
}

var vis = {
  min: 0,
  max: 1,
  palette: [
      'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
      '74A901', '66A000', '529400', '3E8601', '207401', '056201',
      '004C00', '023B01', '012E01', '011D01', '011301'
  ]
};

Map.setCenter(-93.7848, 30.3252, 11);

Map.addLayer(collection.map(NDVI).mean(),vis,'Mean NDVI')
Map.addLayer(collection.map(SAVI).mean(),vis,'Mean SAVI')

Phase2

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_TOA')
      .filterDate('2002-11-01','2002-12-01')

ImageCollectionにて、LANDSAT-7データの2002年11月から12月のデータを取得しています。なお、場所についてはフィルタリングしていないです。

var NDVI = function(image){
  return image.expression('float(b("B4") - b("B3")) / (b("B4") + b("B3"))')
}

var SAVI = function(image){
  return image.expression(
    '(1 + L) * float(nir - red) / (nir + red + L)',
    {
      'nir':image.select('B4'),
      'red':image.select('B3'),
      'L':0.2
    })
}

expressionにてNDVIとSAVIを計算しています。ここではnormalizedDifferece関数は使っていないです。SAVIとは関数が容易されてないため、expressionを使わないといけません。NDVIともSAVIもexpressionとすることで可読性を高めているのだと思われます。

var vis = {
  min: 0,
  max: 1,
  palette: [
      'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
      '74A901', '66A000', '529400', '3E8601', '207401', '056201',
      '004C00', '023B01', '012E01', '011D01', '011301'
  ]
};

ここではmin,max,paletteを設定しています。
書いてて思ったのですが、Paletteの16進数を書くの大変ですよね。
paletteを簡単に設計できるツールがあればすごく便利ですね。

Map.setCenter(-93.7848, 30.3252, 11);

Map.addLayer(collection.map(NDVI).mean(),vis,'Mean NDVI')
Map.addLayer(collection.map(SAVI).mean(),vis,'Mean SAVI')

addLayerにて表示させています。

Filtered Composite

ee.ImageCollection()で複数データを取得し、設定したポリゴン領域のみ表示させます。

Phase1

以下コードです。

var polygon = ee.Geometry.Polygon({
  coords: [[[-109.05, 37.0], [-102.05, 37.0], [-102.05, 41.0], // Colorado
            [-109.05, 41.0], [-111.05, 41.0], [-111.05, 42.0], // Utah
            [-114.05, 42.0], [-114.05, 37.0], [-109.05, 37.0]]],
  geodesic: false
});

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1')
        .filterDate('2000-04-01','2000-07-01')
        .filterBounds(polygon)

var median = collection.median()
var Params = {gain:[1.4,1.4,1.1]}
var result = median.select('B3','B2','B1')
Map.addLayer(result,Params)
Map.setCenter(-110,40,5)

Phase2

var polygon = ee.Geometry.Polygon({
  coords: [[[-109.05, 37.0], [-102.05, 37.0], [-102.05, 41.0], // Colorado
            [-109.05, 41.0], [-111.05, 41.0], [-111.05, 42.0], // Utah
            [-114.05, 42.0], [-114.05, 37.0], [-109.05, 37.0]]],
  geodesic: false
});

ここで領域指定しています。なお、デフォルトでEPSG:4326となっているみたいです。
ポリゴンの複数与えてポリゴンを作成することができます。
coordsは、多角形のポリゴンを定義するためのリストを指定します。
並びは[lon,lat]です。

緯度経度で示されてもわかりにくいので、Google Earthにてピンをつけてみました（すみません。ピンが5つしかつけれず、中途半端になっていますが、なんとなく範囲がわかるかと思います...）

geodesic: falseとあります。地球は楕円形を形状をしていますが、地図の投影法によっては平面になったり、球状になったり、様々あります。いずれにせよ
。地球表面はまっすぐではなく、湾曲しています。geodesicがFalseの場合は平面と仮定し、定義したpolygonをまっすぐ線を引きます。trueの場合は湾曲して線をひきます。
今回の場合は、どちらに設定しても結果は同じ画像が見られます。

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1')
    .filterDate('2000-04-01', '2000-07-01')
    .filterBounds(polygon);

ee.ImageCollectionより、LANDSAT-7のデータを取得しています。期間は2000年4月1日から7月1日の3ヶ月間です。範囲は、Polygonで決定した範囲を指定しているため、全球レベルでのデータは表示されません。

var median = collection.median()
var Params = {gain:[1.4,1.4,1.1]}
var result = median.select('B3','B2','B1')
Map.addLayer(result,Params)
Map.setCenter(-110,40,5)

.selectにてB3、B2、B1を取得し、gainを設定して表示させています。

Linear Fit

DMSPのセンサーにて取得したデータから夜間照明を表示します。

Phase1

以下コードです。

function createTimeBand(img){
  var year = img.date().difference(ee.Date('1999-01-01'),'year')
  return ee.Image(year).float().addBands(img);
}

var collection= ee.ImageCollection('NOAA/DMSP-OLS/CALIBRATED_LIGHTS_V4')
          .select('avg_vis')
          .map(createTimeBand)
var fit = collection.reduce(ee.Reducer.linearFit())

Map.addLayer(ee.Image(collection.select('avg_vis').first()),
      {min:0,max:63},
      'stable lights first asset')

Map.setCenter(30,45,4)
Map.addLayer(fit,{min:0,max:[0.18,20,-0.18],bands:['scale','offset','scale']},
      'stable lights trend')

Phase2

function createTimeBand(img){
  var year = img.date().difference(ee.Date('1999-01-01'),'year')
  return ee.Image(year).float().addBands(img);
}

.date()という関数がでてきています。これは、画像の取得時刻をDateオブジェクトで返す関数です。
その後、Difference関数がでてきています。これは、2つの日付の差を返す関数です。
この場合、imgでのデータの日付と、'1999-01-01'の差を返しています。なお、単位は'year'です。
'year'以外は、'month' 'week', 'day', 'hour', 'minute', or 'second'.があります。

その後、returnにて返していますが、中がee.Image(year).float().addBands(img);とあります。
addBands()の説明は以下の通りです。
「最初の入力からコピーされたすべてのバンドと2番目の入力から選択されたバンドを含むイメージを返します。必要に応じて、最初のイメージのバンドを同じ名前で上書きします。新しい画像は、最初の入力画像からのメタデータとフットプリントを持ちます。」

引数は以下のとおりです。
第一引数；
バンドをコピーする画像
第二引数；
コピーするバンドを含む画像
第三引数；
第二引数で指定した画像が複数バンドを含む場合はそのバンド名を指定（省略した場合は第二引数でした画像のバンドすべてをコピーする）
第四引数；
ここは論理値になります。Trueの場合は第1引数で指定した名前のバンドをすべて上書きします。False場合はリネームします（例えば「foo」という名前の場合は「foo_1」のようにリネームします。

今回の場合は、1996年から1999年までの1年ごとのデータをcollectionしていると思われます。（いや違うかも。誰か指摘お願いします）

var collection= ee.ImageCollection('NOAA/DMSP-OLS/CALIBRATED_LIGHTS_V4')
          .select('avg_vis')
          .map(createTimeBand)

ここではNOAAのDMSPデータを呼び出し、'avg_vis'を指定してcreateTimeBand関数を呼び出しています。
mapは作成したcreateTimeBandのアルゴリズムをcollectionに適用する関数です。

var fit = collection.reduce(ee.Reducer.linearFit())

reduceはImageCollectionや観測時間やバンド、配列といった様々な構造化されたデータを集計するための関数です。例えば最大値、最小値、中央値、標準偏差以外にも、ヒストグラムや線形回帰なども指定できます。
今回の場合、線形回帰をし、夜間光のトレンドをみています。

Map.addLayer(ee.Image(collection.select('avg_vis').first()),
      {min:0,max:63},
      'stable lights first asset')

Map.setCenter(30,45,4)
Map.addLayer(fit,{min:0,max:[0.18,20,-0.18],bands:['scale','offset','scale']},
      'stable lights trend')

下記に、stable lights first assetを示します。

下記にstable lights trendを示します。

今回はサンプルコードにgain値が指定されていたのでいいですが、自分で何か同様のことをやろうとしたときのパラメータ調節が非常に難しいですね。

Simple Cloud Score

ee.ImageCollectionから最も曇りの少ないピクセルを選択し、Landsat8で曇りのない画像をコンポジットします。

Phase1

下記はコードです。

var LC8_BANDS = ['B2','B3','B4','B5','B6','B7','B10']
var STD_NAMES = ['blue','green','red','nir','swir1','swir2','temp']

var cloudScore = function(img){
  var rescale = function(img,exp,thresholds){
    return img.expression(exp,{img: img})
              .subtract(thresholds[0]).divide(thresholds[1] - thresholds[0])
  }

  var score = ee.Image(1.0)
  score = score.min(rescale(img,'img.blue',[0.1,0.3]))

  score = score.min(rescale(img,'img.red + img.green + img.blue',[0.2,0.8]))

  score = score.min(
    rescale(img, 'img.nir + img.swir1 + img.swir2',[0.3,0.8]))

  score = score.min(rescale(img,'img.temp',[300,290]))

  var ndsi = img.normalizedDifference(['green','swir1'])
  return score.min(rescale(ndsi,'img',[0.8,0.6]))
}

var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
                  .filterDate('2017-05-01','2017-07-01')
                  .map(function(img){

                    var score = cloudScore(img.select(LC8_BANDS,STD_NAMES))
                    score = ee.Image(1).subtract(score).select([0],['cloudscore'])
                    return img.addBands(score)
                  })

var vizParams = {bands: ['B4', 'B3', 'B2'], max: 0.4, gamma: 1.6};
Map.setCenter(-120.24487,37.52280,8)
Map.addLayer(collection.qualityMosaic('cloudscore'),vizParams)

Phase2

var LC8_BANDS = ['B2','B3','B4','B5','B6','B7','B10']
var STD_NAMES = ['blue','green','red','nir','swir1','swir2','temp']

センサー名とバンド名を対応付けています。

var cloudScore = function(img) {
  // A helper to apply an expression and linearly rescale the output.
  var rescale = function(img, exp, thresholds) {
    return img.expression(exp, {img: img})
        .subtract(thresholds[0]).divide(thresholds[1] - thresholds[0]);
  };

  // Compute several indicators of cloudyness and take the minimum of them.
  var score = ee.Image(1.0);
  // Clouds are reasonably bright in the blue band.
  score = score.min(rescale(img, 'img.blue', [0.1, 0.3]));

  // Clouds are reasonably bright in all visible bands.
  score = score.min(rescale(img, 'img.red + img.green + img.blue', [0.2, 0.8]));

  // Clouds are reasonably bright in all infrared bands.
  score = score.min(
      rescale(img, 'img.nir + img.swir1 + img.swir2', [0.3, 0.8]));

  // Clouds are reasonably cool in temperature.
  score = score.min(rescale(img, 'img.temp', [300, 290]));

  // However, clouds are not snow.
  var ndsi = img.normalizedDifference(['green', 'swir1']);
  return score.min(rescale(ndsi, 'img', [0.8, 0.6]));
};

クラウドスコアを計算しています。
rescaleは何をしているのかよくわかりませんが、直線に変換？している？だめだよくわからんです。

その後、たくさんのscoreがありますが、各波長に対する雲の応答特性をみてscore値を計算しています。

var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
                  .filterDate('2017-05-01','2017-07-01')
                  .map(function(img){

                    var score = cloudScore(img.select(LC8_BANDS,STD_NAMES))
                    score = ee.Image(1).subtract(score).select([0],['cloudscore'])
                    return img.addBands(score)
                  })

ee.ImageCollectionにてデータを呼び出し、先程作成したcloudScoreを使って雲が最も少ないバンドを作成し、追加しています。

var vizParams = {bands: ['B4', 'B3', 'B2'], max: 0.4, gamma: 1.6};
Map.setCenter(-120.24487,37.52280,8)
Map.addLayer(collection.qualityMosaic('cloudscore'),vizParams)

下記が結果となります。

どうしてこのコードで雲が少ない画像ができあがるのか、いまいち理解しておりません。
またあとで読み直します。

Landsat Simple Composite

Landsatデータを合成します。

Phase1

以下コードです。

var L8 = ee.ImageCollection('LANDSAT/LC08/C01/T1')

var composite = ee.Algorithms.Landsat.simpleComposite({
      collection:L8.filterDate('2015-1-1','2015-7-1'),
      asFloat:true
})

Map.setCenter(-47.6735, -0.6344, 12);
Map.addLayer(composite,{bands:['B6','B5','B4'],max:[0.3,0.4,0.3]})

Phase2

var composite = ee.Algorithms.Landsat.simpleComposite({
      collection:L8.filterDate('2015-1-1','2015-7-1'),
      asFloat:true
})

今回新しくee.Algorithms.Landsat.simpleComposite()がでてきました。
これはすごいです。生データ（RAWデータ）のlandsatデータのCollectionから、雲ができるかぎり入っていないピクセルを採用し、データを作成します。
この一つ前でやったことが、この1行で終わりますね。

Map.setCenter(-47.6735, -0.6344, 12);
Map.addLayer(composite,{bands:['B6','B5','B4'],max:[0.3,0.4,0.3]})

下記に画像を示します。

おわりに

今回は大量にデータを扱うee.ImageCollectionを扱いました。自分のローカルマシンでこのような大量処理を使用とする大量にメモリを使いますし、そもそも衛星データのダウンロードが大変のため、かなり便利だと思います。
個人的には、ee.ImageCollectionにて呼び出したデータを、一枚一枚表示できるようになると嬉しいなと思いました。そもそもそういったことができるかわかりませんが、ee.ImageCollectionにて表示させると時系列関係なく表示されてしまうので、見たい日付が決まっている場合には使いにくいなと思いました。もちろんそういう使い方をする場合はImage.Collectionは使わないほうがいいと思います。

[GEEネタ]Google Earth Engineのサンプルコード写経となんちゃって解説~ImageCollection編~

はじめに

サンプルコードについて

ImageCollection

Clipped Composite

Phase1

Phase2

Expression Map

Phase1

Phase2

Filtered Composite

Phase1

Phase2

Linear Fit

Phase1

Phase2

Simple Cloud Score

Phase1

Phase2

Landsat Simple Composite

Phase1

Phase2

おわりに