クリスマスの聖地を人工衛星から捉える

[Ateam Lifestyle Advent Calendar 2019](https://qiita.com/advent-calendar/2019/ateam-lifestyle) の12日目は株式会社エイチームライフスタイルでWebエンジニアをしている [@water_resistant](https://qiita.com/water_resistant) が担当します。

## はじめに

12月といえばクリスマスがあり、街中光がキラキラして楽しくなる季節ですね。
ただクリスマスは1日しかないので、どうせなら一番キラキラしている場所*(聖地と定義)*で楽しみたいかと思います。
今回はそんな人のためにIT技術を使って**文字通り最高にキラキラ**している場所を探していきます。

## 使ったツール
今回は一番**キラキラ**している場所を探すために、[Google Earth Engine](https://earthengine.google.com/)を使いました。

## GoogleEarthEngineとは
公式ページの紹介をみてみます

> Google Earth Engineは、衛星画像と地理空間データセットの数ペタバイトのカタログと惑星規模の分析機能を組み合わせて、科学者、研究者、開発者が変化を検出し、傾向を地図化し、地球表面の差異を定量化できるようにします
> (Google翻訳済)

つまり衛星画像をデータソースとして様々な分析をできるプラットフォームというわけですね。
今回はこのEarthEngineを**キラキラスポット**探索の為に利用します。

## GooleEarthEngineを使う準備

[SingUpフォーム](https://earthengine.google.com/signup/)にメールアドレスなどの情報と使用用途を記載するだけです。
今回の目的には「明かりをたくさん使っている=電気を沢山使っている場所を調査する」ためと記述。
<img width="767" alt="68747470733a2f2f71696974612d696d6167652d73746f72652e73332e61702d6e6f727468656173742d312e616d617a6f6e6177732e636f6d2f302f363739392f32353739323731362d383033322d366435392d346461342d3263646465636438363461632e706e67.jpg" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/3770c232-e4d4-4ef5-612e-ba7927906ade.jpeg">


1日程度時間が空いてから登録完了のメールがきました。
これを受け取った後でないとAPIを呼び出しても認証エラーになります。↓
<img width="671" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/2eae9590-3752-0151-b445-f3434c2d3e94.png">

実はこれだけで衛星画像へのアクセスや加工が簡単にできるようになります。
**すごいですね。**

今回は触るのが初めてということもあるので、APIではなくIDEからおこなってみます。

IDEはこんな感じです↓
<img width="1792" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/b75f6675-2bdf-c979-d60c-d576f30773b3.png">


## 今回使用する人工衛星データ
IDEを表示できたので、**キラキラスポット**を探すために利用するデータを探します。
[DataCatalog](https://developers.google.com/earth-engine/datasets/catalog)で「light」と検索をしたところ、いくつかデータを発見。
![developers_google_com_earth-engine_datasets_catalog.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/9202088c-ac65-409d-2746-97030af34c60.png)

今回の用途に適してそうなデータ↓
[VIIRS Stray Light Corrected Nighttime Day/Night Band Composites Version 1](https://developers.google.com/earth-engine/datasets/catalog/NOAA_VIIRS_DNB_MONTHLY_V1_VCMSLCFG)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/2fa3bbc4-336b-3069-9a40-0f1d0fb841d8.png)

- 衛星名称
    - [Suomi NPP](https://www.restec.or.jp/satellite/suomi-npp)
- 使用機材
    - Visible/Infrared Imager and Radiometer Suite(VIIRS) / マルチチャンネルイメージャ・放射計

規約についても問題なさそうなのでこれに決定。

> NOAAのデータ、情報、および製品は、配信方法に関係なく、著作権の対象ではなく、その後の一般の使用に対する制限はありません。一度取得すると、それらは合法的に使用できます。前述のデータはパブリックドメインにあり、使用と配布に制限なく提供されています。
> (Google翻訳済)

## データを絞り込むための事前準備

**キラキラスポット**を探すためには、**市区町村(スポット)で区切られたデータ**で絞り込む必要があります。
そのために**行政区域データ**をダウンロードしてEarthEngineに取り込みます。

規約的にも問題なく信用できそうなデータが[政府統計GISデータダウンロード](https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2015&serveyId=A002005212015)にあったので今回はこれを利用。
([いい感じのサイト](https://www.esrij.com/products/japan-shp/)が他にありましたが [規約](https://www.esrij.com/cgi-bin/wp/wp-content/uploads/2015/08/shikuchoson_terms.pdf)的に今回の用途では無理そうなので断念)


あとはデータをダウンロードしてからEarthEngineのIDEからアップロードするだけです。
データが多すぎるのも困るので今回は愛知県のデータのみアップロードします。 (文字コードのみ注意)
<img width="580" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/b55d3a75-6f5e-3c3a-d2c3-eea9eb93d348.png">

これでプログラムから「区画区切りデータ」を利用することができるようになりました。


## 実際にデータを取得する
[公式ドキュメント](https://developers.google.com/earth-engine) などを見ながらコードを書いていく

今回は2014年~2018年の年別データを取得して、光量が多いところや増加傾向のところを**キラキラスポット**とする。

```javascript
// 衛星データのBAND: avg_radから日付を指定してImageCollection取得
var y_2014 = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG').select('avg_rad').filter(ee.Filter.date('2014-12-01T00:00:00+09:00','2014-12-01T23:59:59+09:00')).median();
var y_2015 = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG').select('avg_rad').filter(ee.Filter.date('2015-12-01T00:00:00+09:00','2015-12-01T23:59:59+09:00')).median();
var y_2016 = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG').select('avg_rad').filter(ee.Filter.date('2016-12-01T00:00:00+09:00','2016-12-01T23:59:59+09:00')).median();
var y_2017 = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG').select('avg_rad').filter(ee.Filter.date('2017-12-01T00:00:00+09:00','2017-12-01T23:59:59+09:00')).median();
var y_2018 = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG').select('avg_rad').filter(ee.Filter.date('2018-12-01T00:00:00+09:00','2018-12-01T23:59:59+09:00')).median();

// IDE上のMap表示のオプション
var nighttimeVis = {min: 0.0, max: 60.0, opacity: 0.8}; 
Map.setCenter(136.8809966,35.174377);
Map.setZoom(10);

// IDE上のMapにデータを表示
Map.addLayer(y_2014.clip(aichi_pref_shape), nighttimeVis, '[YEAR] 2014');
Map.addLayer(y_2015.clip(aichi_pref_shape), nighttimeVis, '[YEAR] 2015');
Map.addLayer(y_2016.clip(aichi_pref_shape), nighttimeVis, '[YEAR] 2016');
Map.addLayer(y_2017.clip(aichi_pref_shape), nighttimeVis, '[YEAR] 2017');
Map.addLayer(y_2018.clip(aichi_pref_shape), nighttimeVis, '[YEAR] 2018');

// データ集約方法定義
var reducers_fnc = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.minMax(),
  sharedInputs: true
}).combine({
  reducer2: ee.Reducer.median(),
  sharedInputs: true
});

// データ集約
function reduce_collection(image) {
  return image.reduceRegions({
    collection: aichi_pref_shape,
    reducer: reducers_fnc,
    scale: 500
  })
}

// データ出力タスク定義 (GoogleDrive出力)
function export_table(table, description) {
  return Export.table.toDrive({
    collection: ee.FeatureCollection(table),
    folder: 'EarthEngineExport',
    description: description,
    selectors: (["KEY_CODE", "CITY", "CITY_NAME", "S_NAME", "JINKO", "SETAI", "mean", "max", "min"])
  })
}
// 定義ファイルを参考
// https://www.e-stat.go.jp/gis/statmap-search/data?datatype=2&serveyId=A002005212015&downloadType=1

// データ出力タスク呼び出し
export_table(reduce_collection(y_2014), 'stray_light_2014');
export_table(reduce_collection(y_2015), 'stray_light_2015');
export_table(reduce_collection(y_2016), 'stray_light_2016');
export_table(reduce_collection(y_2017), 'stray_light_2017');
export_table(reduce_collection(y_2018), 'stray_light_2018');
```

上記コードを実行すると、愛知県区画に絞られたデータが結果に出てきます。↓
<img width="1526" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/1c6b47e2-6b64-353f-c018-4c43cac08b62.png">

データ集約も問題なさそうですね。↓
<img width="508" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/291b9573-d20f-b918-87a0-51d396b6adbd.png">

あとはIDEの「Task Run」からデータ出力用の関数を実行します。

GoogleDriveに格納されているのを確認できました。↓
<img width="1421" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/15ff3c7f-2438-2e23-25ff-ae6993f9fedd.png">

IDEはグラフィカルに見ることができるのがいいですが、定量的に測る際は適切ではないので今回はCSV出力してそこでデータを比較することにしました。

## いざデータ加工
出力されたCSVファイルを結合して年別推移を見られるようにしたのが下のキャプチャです。
![1E56FB0F-E421-4D60-A8F8-3E1E1246A6BA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/f829e7ae-5962-3c5e-643c-3658d64d66bc.png)

[CAGR](https://mba.globis.ac.jp/about_mba/glossary/detail-11621.html)でソートしてみたところ、「岡崎市 宮石町」が2016年から4倍跳ね上がっていたので調べてみます。
<img width="1549" alt="スクリーンショット 2019-12-09 1.02.35.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/1b63bd99-f732-edf3-d727-a192842ecf92.png">

該当住所をGoogleMapで調べたところ、エリア内にSA(サービスエリア)を見つけました。↓
<img width="1175" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/45543ac9-cb61-14b8-d799-e82648aa60f2.png">

[NEXCO中日本、2016年2月13日15時にオープンする新東名「岡崎SA」](https://travel.watch.impress.co.jp/docs/news/742997.html%0A) 
なるほど面白いですね。

ただ[CAGR](https://mba.globis.ac.jp/about_mba/glossary/detail-11621.html)ランキング上位を見ても、私が求めている**キラキラ**とは少し違うので、別の観点「過去5年間の平均値」で見ていきます。

![F8B2A403-67C4-4EE9-AF56-554642B88955.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/f4bf568a-4fa5-29cc-d8b7-209d93f92a9d.png)

上位5件はこうなりました。

1. 中区 錦３丁目
2. 中区 新栄町１丁目
3. 中区 栄４丁目
4. 一宮市 丹陽町三ツ井字平山
5. 中区 新栄町２丁目

|中区 錦３丁目|中区 新栄町１丁目|中区 栄４丁目|一宮市 丹陽町三ツ井字平山|中区	新栄町２丁目|
|:--:|:--:|:--:|:--:|:--:|
| <img width="616" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/5ec8c08c-494f-6d4b-2f30-e89119134ed0.png">|<img width="420" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/798dda25-2b8b-05b1-4701-698fc196d3b8.png">|<img width="468" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/7eae6c98-dcad-42f9-a896-0f3cff18b70e.png">|<img width="910" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/cac9822d-7538-8549-e2c8-6f2d92c0b34b.png">|<img width="657" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/c089a4c8-4ec3-b2cd-22f6-073f1aa3f9fd.png">|


「栄」や「錦」が出ることは予想してましたが、、 4位は完全に予想外です。何か変な建物が見えますが...。
(実は一宮は4位だけでなく、それ以降もそれなりにでてきます。)

ただひょっとしたら**キラキラスポット(聖地)**の答えがここにあるかもしれないと思ったので少し探ります。

衛星の分解能が400m-800mで、区間分割がうまく行かずに明るくなっていると思われるので
該当の住所をEarthEngine側で拡大してみましょう。

<img width="1243" alt="pic" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6799/0d9bd65b-e873-7d3c-fa30-28ba5552da41.png">


**IC(インターチェンジ)近く**が無茶苦茶明るい様子。
確かに豪華絢爛なお城が多いですもんね。
知りたくなかった。

## まとめ

GoogleEarthEngineは、水温や気温を始めとして凄く沢山のデータがあります。
例えば下のようなデータについては図書館に行かずともPCの前に居ながら出すことができます。

- 森林がどれだけ減っているのか
- 気温/海温はどれだけ変動しているのか

使い方は簡単ながらも、非常に多くを知ることができるので、自分の視点で地球を一度見てみるのも面白いかと思います。
(本来はその用途のツールです)

あとクリスマスに向けての**キラキラスポット(聖地)調査**は全くもって役に立たなかった。

---

 [Ateam Lifestyle Advent Calendar 2019](https://qiita.com/advent-calendar/2019/ateam-lifestyle)  の 13日目は、 [@ponco](https://qiita.com/ponco)がお送りします!
私のネタエントリとは違う、しっかりした記事を書いてきてくれることでしょう!

そんな”挑戦”を大事にするエイチームグループでは、一緒に働けるチャレンジ精神旺盛な仲間を募集しています。興味を持たれた方はぜひ[エイチームグループ採用サイト](https://www.a-tm.co.jp/recruit/)を御覧ください。


## 今回のコード
今回作成したデータやコードへのリンクを貼っておきます。
(EarthEngineはアカウント登録してないと見れないです。)

- [GoogleEarthEngine Code](https://code.earthengine.google.com/cac84e15e6fdcfac4b8ddade4fc1e095)
- [stray_light_combine - Google スプレッドシート](https://docs.google.com/spreadsheets/d/1FwxLYty8-evnDZTARgfE_ue8kGGNHGpArA8OMcS0Z98/edit?usp=sharing%0A)


## 参考・データ元

- [衛星画像解析が変わる!? 「Google Earth Engine」の何がすごいのか | 宙畑](https://sorabatake.jp/337/)
- [統計GISデータダウンロード | 政府統計の総合窓口](https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2015&serveyId=A002005212015)
- [Google Japan Blog: Google Earth Engineが可能にする新しい Google Earth タイムラプス](https://japan.googleblog.com/2019/04/google-earth-engine-google-earth.html)
- [GEEネタGoogle Earth Engineのサンプルコード写経となんちゃって解説~Image編~ - Qiita](https://qiita.com/iwasaki_kenichi/items/0f8087b870688916525e)