はじめに
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つを自分に課したいと思います。
■Phase1:まずは普通に写経(コピペは厳禁)
■Phase2:コードの解説
※Image編では、
■Phase3:写経したコードを元に少し改良する(余裕があったら)
を入れていましたが、外します。ちょっと時間がかかる...
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を作成し、似たようなことをしてみたいと思いました。