QGIS上で道路周辺の夜間光の変化を見る

# 1. はじめに
今回の分析対象はネパールのシンズリ道路です。
シンズリ道路については、[JICAプロジェクト概要ページ](https://www.jica.go.jp/oda/project/1860370/index.h!)をご覧ください。

この道路の建設によって、対象地域の「経済活性化」が期待されていました。
そこで、衛星データ（夜間光）を活用して、道路による「経済活性化」を可視化してみたいと思います。

今回使用するツールはGoogle Earth Engine（以下、「GEE」）とQGISです。
GEE上で衛星データのインポート・成形・ダウンロードを行い、QGISでそのダウンロードしたデータを集計・可視化します。

最終的なアウトプットは下図になります。
このように図示することによって、どこでどのような変化があったのかがわかりますね。
![ntl_change_nepal.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/5550b2b3-f034-b244-ac0c-d1f80c938db4.png)

さらに、シンズリ道路周辺に注目すると、道路沿いの町（特にドリケルやシンズリマディ）において夜間光量が増加しており、事業が町の「経済活性化」に貢献していることが示唆されます。
![ntl_change_road.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/b6e73326-2c78-90ae-7557-4059c55f4e06.png)


以下、上記の図を作成するための手順を紹介します。

# 2. 使用データ
今回の分析には、以下の衛星データを使います。
- 夜間光：[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)
    - データ期間：2014～2023（2024年4月時点で2023年11月までのデータが利用可能）
    - 解像度：約500m

また、衛星データに加えて、以下の地理データも使います。
- 行政区：[Nepal - Subnational Administrative Boundaries](https://data.humdata.org/dataset/cod-ab-npl?)
    - adm1：州レベル、adm2：県レベル、adm3：区レベル
    - データ最終更新日：2024年3月<br><br>

- 道路と町：無料公開データがないため、自作で対応します。
    - QGIS上でGoogle Mapのレイヤーを重ねて目視で道路と町を特定し、データを自ら作成します。<br>
    - 作成方法については、[こちら（過去の類似記事）](https://qiita.com/ishimotojuri/items/fcdb91cf8f4601f38d52)をご覧ください。

# 3. 分析手順
## 3.1.  衛星データの収集
Google Earth Engineを使って、夜間光データを収集します。

<INS>まず、分析領域を設定します</INS>。
国境データをGEEにインポートし、ネパールのみを選択します。
<details><summary>インポートのコードはこちら</summary><div>

```js:sample.js
var dataset = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');//境界データを取得する
var nepal = dataset.filter(ee.Filter.eq('country_na', 'Nepal'));//ネパールのデータを取得する
// ('シェープファイル上の変数名', 'データ名')

Map.addLayer(nepal, {
  fillColor: 'b5ffb4',
  color: '00909F',
  width: 1.0,
},
'青:ネパール');
// {}：地図に表示される方法（色など）。デフォルトは黒

```
</div></details>

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/344e2689-810c-ab2f-2b6e-6cb8d6c2f7c4.png)


<INS>次に、事業実施前後の夜間光データをインポートします</INS>。
各条件を次のように設定し、GEEに夜間光データをインポートします。

①事業実施前（2014年）
データ期間：2014-04-01, 2024-12-31
データ種類：avg_rad
分析領域：nepal（上記で設定した範囲）
統計量：mean（平均値を算出）

②事業実施後（2023年）
データ期間：2023-01-01, 2023-12-31
データ種類以下、同上

<details><summary>インポートのコードはこちら</summary><div>

```js:sample.js
var pre_data = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG") 
             .filterDate('2014-01-01', '2014-12-31')
             .select('avg_rad')
             .filterBounds(nepal)
             .mean();

var post_data = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG") 
             .filterDate('2023-01-01', '2023-12-31')
             .select('avg_rad')
             .filterBounds(nepal)
             .mean();

```
</div></details>

<INS>最後に、夜間光の画像をエクスポートします</INS>。
ネパールの国境に沿った画像にして（切り抜き）、エクスポートします。

<details><summary>エクスポートのコードはこちら</summary><div>

```js:sample.js
// pre data
var pre_data_clip = pre_data.clip(nepal);

Export.image.toDrive({
  image: pre_data_clip,
  description: 'ntl_nepal_2014',
  folder: 'Nepal', // set a folder of your google drive
  region: nepal, 
  scale: 500, //resolusion of the image
  crs: 'EPSG:4326'
});

// post data
var post_data_clip = post_data.clip(nepal);

Export.image.toDrive({
  image: post_data_clip,
  description: 'ntl_nepal_2023',
  folder: 'Nepal', // set a folder of your google drive
  region: nepal, 
  scale: 500, //resolusion of the image
  crs: 'EPSG:4326'
});
```
</div></details>

上記のコードを実行すると、右側の「Tasks」画面で各ファイル名が表示されます。
ファイル名の「Run」ボタンをクリックすると、画面が表示されます。

その画面の設定条件を確認し、下の「Run」ボタンをクリックするとエクスポートが開始されます。

ダウンロードしたデータはGoogleドライブの任意のフォルダに保存されます。

## 3.2. データ作成
<INS>分析領域を作成します</INS>。
上記でダウンロードした行政区、作成した道路と町の地理データをQGISのメイン画面にドラッグ＆ドロップして表示させます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/58bbe25f-fc5a-79ec-a34e-60c8b77c8c1f.png)

<INS>夜間光データを追加します。</INS>
ダウンロードした二枚の画像（事前と事後）をQGISのメイン画面にドラッグ＆ドロップします。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/f02f409e-a2fe-9006-86ba-c693400b8898.png)

<INS>夜間光量の変化を計算します。</INS>
「ラスタ」→「ラスタ計算機」を選択します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/20b2b56c-a38c-c4ff-0ffd-49beaa7d0f0a.png)

以下の通り設定し、「OK」をクリックします。
式："ntl_nepal_2023@1" - "ntl_nepal_2014@1"
バンド、演算子をそれぞれダブルクリックすると、式に反映されます。
出力レイヤ："diff_ntl_nepal"という名前で任意のフォルダに保存
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/8699851e-51b6-baad-63fd-50e91d103c56.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/2aab71ad-a346-383f-8513-52e78da26e37.png)

<INS>算出した夜間光量の変化を可視化します。</INS>
新しく追加された「diff_ntl_nepal」レイヤを右クリック→「プロパティ」→「シンポロジ」を選択します。

以下の通り設定し、「適用」→「OK」をクリックします。
+ レンダリングタイプ：「単バンド疑似カラー」
+ カラーランプ：Spectral→その上にある「カラーランプを反転」を選択
+ 分類モード：連続的
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/4b910554-7e41-6c74-5dbe-6e8514738281.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/c385b396-407a-c3c0-89a2-3d0d74caea85.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/28f616e9-5bae-7316-7b2e-b3d20ac215ef.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/b1787753-15c1-515d-1aee-84bff98d7fc7.png)

完了後の画像は以下の通りです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/712d4723-285c-5a39-87d6-a874af779bcf.png)

なお、最小/最大範囲について、デフォルトは「範囲全体のすべての値を使用する」となっています。しかし、ごく少数の非常に小さい値や非常に大きい値があると、ラスタのレンダリングに悪影響を及ぼすことがあります。 外れ値があると、画像を可視化した際に、データの全体的な傾向やパターンが正しく表現されなくなります。

「ヒストグラム」でデータの分布を見たところ、ほとんどのピクセルが0-1の間にあり、最小値や最大値はほとんどないことが確認できます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/b7dc7503-95a5-be08-c4ef-7f5484c7364d.png)

そのため、その影響を取り除くために、「最大／最小範囲」を以下の通り変更しました。

最小／最大範囲：平均 +/- 標準偏差 x2.0
→全体の約95％の範囲内の値のみを考慮したカラーテーブルを作成します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/494c9119-2b57-29a5-c3d0-6d24ba3f7725.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/8075b7b6-9250-6d84-1348-39242681332c.png)

## 3.3. 地図の作成
<INS>地図のレイアウト設定を行います。</INS>
「プロジェクト」→「新規印刷レイアウト」を選択し、印刷レイアウトのタイトル（好きな名前で構いません）を入力すると、地図の作成画面が現れます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/9560015d-fe42-7952-f68a-559339f5ffaf.png)

＜地図レイアウト＞
左のツールバーにある「新しい地図をレイアウトに追加」をクリックして、キャンバス上で適当にドラッグし、地図を追加します。
さらに、以下の通り設定し、緯度経度を追加します。
+ 右下ボックスの「アイテムプロパティ」→「フレーム」にチェック
+ 右下ボックスの「アイテムプロパティ」→「グリッド」→「+」（新規グリッドの追加）→作成したグリッド1を選択した状態で「グリッドの修正」を選択
    + グリッドタイプ：フレームと注記のみ
    + 座標参照系：「EPSG：4326」
    + 間隔： XとYどちらも5.0
    + フレームスタイル：外側に目盛り
        + 右側と上側のチェックを外す
    + 「座標を描画」にチェック
        + 左：水平
        + 右：無効
        + 上部：無効
        + 下部：水平
        + フォント：Arial, サイズ9
        + 座標精度 → 0

＜凡例＞
左のツールバーにある「「新しい凡例をレイアウトに追加」をクリックして、キャンバス上で適当にドラッグし、凡例を追加します。さらに、以下の通り設定します。
+ レイヤ上の凡例をクリック→右側の「凡例アイテム」→「リンク先地図の内側にあるアイテムだけ表示する」
+ 「自動更新」のチェックを外す→必要な凡例以外を削除（―）
+ diff_ntl_nepalのカラーパレットをダブルクリック→値を小数点第二位切捨てに修正
+ 「フレーム」にチェック

＜その他＞
左のツールバーにある「新しいスケールバーをレイアウトに追加」、「新しい方位記号をレイアウトに追加」をクリックして、それぞれキャンパスに適当にドラッグし、追加します。方位記号については、▼SVGパラメータ→塗りつぶし色：黒を選択します。

## 3.4. 地図のエクスポート
地図が完成したら、「レイアウト」→「画像としてエクスポート」を選択します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/42f0a5c8-71b4-058b-f050-dce61a1e199d.png)

# まとめ
+ この記事では、衛星・地理データを活用して、事業による「経済活性化」を可視化してみました。
+ 事業実施前と比較すると、事業実施後は、シンズリ道路沿いの町（特にドリケルやシンズリマディ）において夜間光量が増加しており、道路が町の「経済活性化」に貢献していることが示唆されます。