はじめに
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ファイルの中で異なる領域と定義されているものをランダムな色でペイントさせている感じでしょうか。
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
})
その後、表示させます。
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です。許容誤差の値が大きくなると、バッファのマージンの重なりがへります。つまり、バッファエリアの形状が変わります。
エリアマーケティングや商圏分析などで使えるのではないでしょうか？
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);
下記に結果を示します。
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を作成します。
以下のような画像が表示できます。
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)
結果を表示します。下記です。
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'});
最後に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')
最後に表示します。下記の通りです。
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万を超えるラインフィーチャが含まれているデータです。
一部ですが、ピンク色の部分が道路になります。
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')
ここで表示しています。
おわりに
FeatureCollectionが終わりました。Google Earth Engineは衛星データ処理だけではなく、膨大な衛星データも使えるwebGISにだなと強く感じました。
今回はGoogle Earth Engine側が用意しているFusion Tableでした。それですと日本国内のデータは扱えないので、日本国内のFusion Tableを作成し、似たようなことをしてみたいと思いました。