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

# はじめに

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


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コースあります。

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

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

# Feature Collection
Feature Collectionは、三番目のコースです。
今までは、`ee.Image`や`ee.ImageCollection`といった画像データを扱ってきました。`Feature`は画像といったラスタデータはなく、ベクタデータとなります。

##From Fusion Table
Google Earth EngineはGoogle Fusion Tableのデータを扱うことができるみたいです。（Google Fusion Tableについてはこちらを参照してください：http://blog.thetheorier.com/entry/fusion-tables）

また、Google Earth Engineで使用するベクターデータセットは、Fusion TableとEarth Engine Catalogがあります。
この2つについては下記を参照してください。
https://developers.google.com/earth-engine/vector_datasets

###Phase1
下記コードです。

```
var fc = ee.FeatureCollection('ft:1Ec8IWsP8asxN-ywSqgXWMuBaxI6pPaeh6hC64lA')
      .filter(ee.Filter.stringContains('ECO_NAME','desert'))

var image1 = ee.Image().byte().paint(fc,'ECO_NUM')

print(fc)

Map.setCenter(-109.687, 30.524, 5);
Map.addLayer(image1.randomVisualizer())


Export.table.toDrive({
  collection: fc,
  description: 'fc_table',
  fileFormat: 'CSV'
});
```

###Phase2

```
var fc = ee.FeatureCollection('ft:1Ec8IWsP8asxN-ywSqgXWMuBaxI6pPaeh6hC64lA')
      .filter(ee.Filter.stringContains('ECO_NAME','desert'))
```
ここでFusion Tableを読み込んでいます。さらに`.filter`で`ee.Filter.stringContains`をしています。
'ECO_NAME'内にある'desert'を読み込むということをしていると思われます。しかし、今回のFusion Tableは用意されているデータであり、どんなデータかわかりません。なので、データをexportしたいと思います。
それが

```
Export.table.toDrive({
  collection: fc,
  description: 'fc_table',
  fileFormat: 'CSV'
});
```
です。
4MBのCSVファイルであり、カラム数もインデックス数もそれなりに多かったので、内容については割愛しますが、ECO_NAMEとdesertという文字列はありました。



```
Map.setCenter(-109.687, 30.524, 5);
Map.addLayer(image1.randomVisualizer())
```
ここで画面に表示させていますが、`randomVisualizer()`というものがでてきました。
これはピクセルの各一意の値にランダムないろを割り当てる関数です。
今回は、衛星画像は読み込んでおらずCSVファイルのみのため、CSVファイルの中で異なる領域と定義されているものをランダムな色でペイントさせている感じでしょうか。
<img width="282" alt="img1.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/103b5ef1-423d-c1f2-395a-53f4e4eb02fa.png">


##From Polygons
ポリゴンから`FeatureCollection`を作成し、画面上に表示します。

###Phase1
下記コードです。

```
var fc = ee.FeatureCollection([
    ee.Feature(
      ee.Geometry.Polygon({
        coords: [[-109.05, 41], [-109.05, 37], [-102.05, 37], [-102.05, 41]],
        geodesic:false,
        maxError:1000
      }),{name:'Colorado',fill:1}),
    ee.Feature(
      ee.Geometry.Polygon({
        coords:[
          [-114.05, 37.0], [-109.05, 37.0], [-109.05, 41.0],
          [-111.05, 41.0], [-111.05, 42.0], [-114.05, 42.0]
        ],
        geodesic:false,
        maxError:1000
        }),{name:'Utah',fill:2})
])

var image = ee.Image().toByte()
    .paint(fc,'fill')
    .paint(fc,3,5)
    
Map.addLayer(image,{
  palette: ['000000', 'FF0000', '00FF00', '0000FF'],
  max:3,
  opacity:0.5
  
})

Map.setCenter(-107, 41, 6);
```

###Phase2
```
var fc = ee.FeatureCollection([
    ee.Feature(
      ee.Geometry.Polygon({
        coords: [[-109.05, 41], [-109.05, 37], [-102.05, 37], [-102.05, 41]],
        geodesic:false,
        maxError:1000
      }),{name:'Colorado',fill:1}),
    ee.Feature(
      ee.Geometry.Polygon({
        coords:[
          [-114.05, 37.0], [-109.05, 37.0], [-109.05, 41.0],
          [-111.05, 41.0], [-111.05, 42.0], [-114.05, 42.0]
        ],
        geodesic:false,
        maxError:1000
        }),{name:'Utah',fill:2})
])
```
`ImageCollection`は`Image`を複数スタックされているデータのことでしたが、`FeatureCollection`も同様で、`Feature`を複数スタックしているデータのようです。
`var fc = ee.FeatureCollection`のなかに

```
 ee.Feature(
      ee.Geometry.Polygon({
        coords: [[-109.05, 41], [-109.05, 37], [-102.05, 37], [-102.05, 41]],
        geodesic:false,
        maxError:1000
      }),{name:'Colorado',fill:1}),
```
これと

```
ee.Feature(
      ee.Geometry.Polygon({
        coords:[
          [-114.05, 37.0], [-109.05, 37.0], [-109.05, 41.0],
          [-111.05, 41.0], [-111.05, 42.0], [-114.05, 42.0]
        ],
        geodesic:false,
        maxError:1000
        }),{name:'Utah',fill:2})
```
これがありますね。つまり`ee.Feature`があります。

上段の`ee.Feature`はコロラド州を示しており、下段の`ee.Feature`はユタ州を示しています。

```
var image = ee.Image().toByte()
    .paint(fc,'fill')
    .paint(fc,3,5)
```
その後、'fill'にて塗りつぶし、色を3、線幅を5として定義します。

```
Map.addLayer(image,{
  palette: ['000000', 'FF0000', '00FF00', '0000FF'],
  max:3,
  opacity:0.5
})
```
その後、表示させます。
<img width="282" alt="img2.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/75bc9192-db9e-f1cd-9e57-5f8607df672c.png">


##Buffer
`ee.Feature`のbufferのサンプルです。
サンフランシスコの鉄道の駅から2km以内の地域を表示します。

###Phase1
下記コードです。

```
var bartStations = ee.FeatureCollection('ft:1xCCZkVn8DIkB7i7RVkvsYWxAxsdsQZ6SbD9PCXw');

var buffered = bartStations.map(function(f){
  return f.buffer(2000,100)
})

Map.addLayer(buffered,{colo:'800080'})
Map.setCenter(-122.4,37.7,11)
```

###Phase2

```
var buffered = bartStations.map(function(f){
  return f.buffer(2000,100)
})
```
バッファとは、いろいろな意味があるかと思います。（電子回路でいうと緩衝増幅器であったり、計算機でいうとメモリ管理であったり）
Google Earth Engine（というよりもGISでは）のバッファは特定の距離の範囲を作成することです。（参考：https://www.pasco.co.jp/recommend/cook/cook062/）

この場合、半径が2000mで許容誤差が100です。許容誤差の値が大きくなると、バッファのマージンの重なりがへります。つまり、バッファエリアの形状が変わります。
<img width="281" alt="img3.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/6702a6c0-6a92-e212-57ca-81b15d00628d.png">

エリアマーケティングや商圏分析などで使えるのではないでしょうか？

##Computed Area Filter
面積を計算します。面積が3平方km未満の米国の都市を探します。

###Phase1
下記コードです。

```
var counties = ee.FeatureCollection('ft:1pjtcfSKIbYbj4wRcBjc0Bb6NB-sQRI-L2nIzHiU');

var countiesWithArea = counties.map(function (f){
  
  var areaHa = f.area().divide(100 * 100)
  
  return f.set({area:areaHa})
  
})

var smallCounties = countiesWithArea.filter(ee.Filter.lt('area',3e5))

Map.addLayer(smallCounties,{color:'900000'})
Map.setCenter(-119.7,38.26,7)
```
###Phase2

```
var counties = ee.FeatureCollection('ft:1pjtcfSKIbYbj4wRcBjc0Bb6NB-sQRI-L2nIzHiU');
```
Fusion Tableを読み込みます。

```
var countiesWithArea = counties.map(function (f){
  var areaHa = f.area().divide(100 * 100)
  return f.set({area:areaHa})

})
```
それぞれの都市の面積計算をする関数です。
面積をヘクタールに変換しています。

`return f.set({area:areaHa})`
これは、プロパティにareaという新しい属性情報を追加しています。中身は、前段処理した、面積（ヘクタール）です。

```
ar smallCounties = countiesWithArea.filter(ee.Filter.lt('area', 3e5));
```
ここでは小さな都市のみを取得しています。

```
Map.addLayer(smallCounties, {color: '900000'});

Map.setCenter(-119.7, 38.26, 7);
```
下記に結果を示します。
<img width="284" alt="img4.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/541f4a05-21da-66f4-b95e-0a46e5be8802.png">


##Distance
`ee.FeatureCollection`内で、もっとも近い`Feature`までの距離を計算します。

###Phase1
下記コードです。

```
var fc = ee.FeatureCollection([
  ee.Geometry.Point(-72.94411, 41.32902),
  ee.Geometry.Point(-72.94411, 41.33402),
  ee.Geometry.Point(-72.94411, 41.33902),
  
  ee.Geometry.LineString(
      -72.93411, 41.30902, -72.93411, 41.31902, -72.94411, 41.31902)
  ])
  
var distance = fc.distance(1000,100)
var vis_Params = {min:0,max:1000,palette:['yellow','red']}
Map.setCenter(-72.94, 41.32, 13);
Map.addLayer(distance,vis_Params)
Map.addLayer(fc)
```


###Phase2

```
var fc = ee.FeatureCollection([
  ee.Geometry.Point(-72.94411, 41.32902),
  ee.Geometry.Point(-72.94411, 41.33402),
  ee.Geometry.Point(-72.94411, 41.33902),
  
  ee.Geometry.LineString(
      -72.93411, 41.30902, -72.93411, 41.31902, -72.94411, 41.31902)
  ])
```
`ee.Geometry.point`と`ee.Geometry.LineString`のリストから`ee.FeatureCollection`を作成します。
以下のような画像が表示できます。
<img width="281" alt="img5.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/3d3ad1fa-691d-08bb-88b4-55bc41f99431.png">

```
var distance = fc.distance(1000, 100);
Map.setCenter(-72.94, 41.32, 13);
```
その後、`fc.distance`をしています。
ピクセルの中心から`ee.FeatureCollection`内の任意の`Point`また`LineString`の最も近い部分までの距離（メートル単位）であるイメージを作成します。

第一引数は、各ピクセルからの最大距離（メートル単位）です。この距離内に`Point`や
`LineString`がない場合、ピクセルはマスクされます。

第二引数は、最大投影誤差です。再投影が必要な場合にのみ使用されます。EPSG等を変更する再投影（reprojection）のとき、微妙にずれるときがあります。そのときに使われるのだと思いますが、詳しくは理解できず...

```
var vis_Params = {min:0,max:1000,palette:['yellow','red']}
Map.setCenter(-72.94, 41.32, 13);
Map.addLayer(distance,vis_Params)
Map.addLayer(fc)
```
結果を表示します。下記です。
<img width="282" alt="img6.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/8541e4b8-0cc2-eb30-f2ed-e8014209f38a.png">



##Join
BARTの駅から2km以内のサンフランシスコにある公園を表示します。

###Phase1
下記コードです。

```
var bart = ee.FeatureCollection('ft:1xCCZkVn8DIkB7i7RVkvsYWxAxsdsQZ6SbD9PCXw');
var parks = ee.FeatureCollection('ft:10KC6VfBWMUvNcuxU7mbSEg__F_4UVe9uDkCldBw');

var joinFilter = ee.Filter.withinDistance({
  distance: 2000,
  leftField:'.geo',
  rightField:'.geo'
})

var closeParks = ee.Join.simple().apply({
  primary:parks,
  secondary:bart,
  condition:joinFilter
})

var bufferedBart = bart.map(function(f){
    return f.buffer(2000,100)})

Map.setCenter(-122.45, 37.75, 13);
Map.addLayer(bufferedBart, {color: 'b0b0b0'});
Map.addLayer(closeParks, {color: '008000'});
```

###Phase2
```
var bart = ee.FeatureCollection('ft:1xCCZkVn8DIkB7i7RVkvsYWxAxsdsQZ6SbD9PCXw');
var parks = ee.FeatureCollection('ft:10KC6VfBWMUvNcuxU7mbSEg__F_4UVe9uDkCldBw');
```
ここで読み込んでいるのは、駅と公園のデータです。
※あくまでもサンプルコードだからいいですが、日本でこういうのをしたい場合は自分でFusion Tableを作らなければいけないのは億劫ですね。すでに用意してあるかもしれませんが...

```
var joinFilter = ee.Filter.withinDistance({
  distance: 2000,
  leftField: '.geo',
  rightField: '.geo'
});
```
ここでは後段処理で使用するフィルタを自作しています。
`leftField`が`rightField`から指定された距離（`distance`のことです）内にある場合のみ通過するバイナリフィルタを作成しています。

```
var closeParks = ee.Join.simple().apply({
  primary:parks,
  secondary:bart,
  condition:joinFilter
})
```
前段で作成した`joinFilter`を適用します。
`primary`は`leftField`に対応し
'secondary'は`rightField`に対応します。
つまり、`bart`から2km以内にある`parks`のみを抽出しています。

```
var bufferedBart = bart.map(function(f){
    return f.buffer(2000,100)})
```
`buffer`を使って周囲2kmをモニタに表示するように`buffereBart`を作成しています。

```
Map.addLayer(bufferedBart, {color: 'b0b0b0'});
Map.addLayer(closeParks, {color: '008000'});
```
<img width="283" alt="img7.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/20331653-355d-e0d9-effc-20e8171af3bf.png">

最後に`Map.addLayer()`にてモニタに表示させます。

結果を見ると、完全に包含されている公園を表示させるわけではないみたいですね
公園の一部でも周囲2kmに触れていれば表示するみたいですね。

##Reduce To Image
`ee.FeatureCollection.reduceToImage()`の使い方を学びます。
`ee.FeatureCollection.reduceToImage()`は、交差するすべての`Feature`の選択されたプロパティに`reducer`を適用することで`FeatureCollection`から画像を作成します。

###Phase1
下記コードです。

```
var fc = new ee.FeatureCollection([
  ee.Feature(
    ee.Geometry.Rectangle(
      -122.4550, 37.8035,
      -122.4781, 37.7935),
      {'value':0}),
  ee.Feature(
    ee.Geometry.Polygon([
      [-122.4427, 37.8027],
      [-122.4587, 37.7987],
      [-122.4440, 37.7934]]),
    {'value': 1})
  ]);

var image_reduced = fc.reduceToImage(['value'],'mean')
var vis_Params = {min:0,max:1,palette:['008800','00FF00']}
Map.setCenter(-122.4561, 37.7983, 14);
Map.addLayer(image_reduced,vis_Params,'res')
```
###Phase2
```
var fc = new ee.FeatureCollection([
  ee.Feature(
    ee.Geometry.Rectangle(
      -122.4550, 37.8035,
      -122.4781, 37.7935),
      {'value':0}),
  ee.Feature(
    ee.Geometry.Polygon([
      [-122.4427, 37.8027],
      [-122.4587, 37.7987],
      [-122.4440, 37.7934]]),
    {'value': 1})
  ]);
```
`ee.FeatureCollection`を作成します。例により、`ee.Geometry.Rectangle`と`ee.Geometry.Polygon`を内部に作成します。また、RectangleとPolygonを区別するために'value'の値をそれぞれ異なる値にします。

```
var image_reduced = fc.reduceToImage(['value'],'mean')
```
`fc`を画像にします。つまり`ee.FeatureCollection`を画像にしています。
画像にした際のピクセル値は、`Feature`が交差する'value'の平均値となります。
今回の場合、0.5となります。**（たぶん）**

```
var vis_Params = {min:0,max:1,palette:['008800','00FF00']}
Map.setCenter(-122.4561, 37.7983, 14);
Map.addLayer(image_reduced,vis_Params,'res')
```
最後に表示します。下記の通りです。
<img width="283" alt="img8.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/d64039a6-db69-4871-b336-79c8dc471fd3.png">


##From Earth Engine Asset
今まで何度かGoogle Fusion Tableを使った処理をしてきましたが、今回はEarth Engine Tableというものを使い、`ee.FeatureCollection`を作成します。

###Phase1
下記コードです。

```
var roads = ee.FeatureCollection('TIGER/2016/Roads');
var interstates = roads.filter(ee.Filter.eq('rttyp','I'))
var surfaceRoads = roads.filter(ee.Filter.eq('rttyp','M'))
Map.addLayer(surfaceRoads,{color:'gray'},'surface roads')
Map.addLayer(interstates,{color:'red'},'intestates')
```

###Phase2
```
var roads = ee.FeatureCollection('TIGER/2016/Roads');
```
ここでEarth Engine Catalogの道路データを読み込んでいます。
なお、今回読み込んでいるデータは米国国税調査局が調査した2016年の道路セグメントです。
米国、コロンビア特別区、プエルトリコ、及び各離島をカバーする1900万を超えるラインフィーチャが含まれているデータです。
<img width="281" alt="img9.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/2d0d9921-fbee-a307-3ac2-f22262b61763.png">
一部ですが、ピンク色の部分が道路になります。

```
var interstates = roads.filter(ee.Filter.eq('rttyp','I'))
```
ここは高速道路のみを抽出しています。

```
var surfaceRoads = roads.filter(ee.Filter.eq('rttyp','M'))
```
ここは地上道路のみを抽出しています。

※なぜこの引数で高速道路なのか、地上道路なのか、という疑問がでてきます。
https://code.earthengine.google.com/dataset/TIGER/2016/Roads
上記サイトを確認してください。
Table SchemaのNameに'rttyp'がありますが、これは道路のタイプを示しています。

https://www.census.gov/geo/reference/rttyp.html
上記サイトを確認すると、道路タイプ別コードがあります。これでなぜ高速道路なのか、地上道路なのかが判別できます。

```
Map.addLayer(surfaceRoads,{color:'gray'},'surface roads')
Map.addLayer(interstates,{color:'red'},'intestates')
```
ここで表示しています。
<img width="282" alt="img10.PNG" src="https://qiita-image-store.s3.amazonaws.com/0/284911/9c1b6d53-9587-8366-6661-6dcb1a64bcc1.png">

# おわりに
`FeatureCollection`が終わりました。Google Earth Engineは衛星データ処理だけではなく、膨大な衛星データも使えるwebGISにだなと強く感じました。
今回はGoogle Earth Engine側が用意しているFusion Tableでした。それですと日本国内のデータは扱えないので、日本国内のFusion Tableを作成し、似たようなことをしてみたいと思いました。