夜間光データで電化率を計算してみた

# はじめに
夜間光データを用いてサブサハラアフリカの電化率を計算した[論文（Falchetta. et al., 2019）](https://www.nature.com/articles/s41597-019-0122-6)を見つけました。これを参照し、今回はイラクを対象に電化率を計算したいと思います。

イラク北部のクルド地域（ドホーク県、エルビル県、スレイマニア県）では電力供給が電力需要に追いついていない状態です。2007年時点で、イラク北部のドホーク県では日平均10～11時間、エルビル県及びスレイマニア県では日平均5時間～6時間程度しか電力供給のない状態でした。

またイラク南部では、2010年及び2022年に8月の気温50度の暑さの中、多くの住民が停電に見舞われたことに抗議し、デモが行われたとの報道がありました。

このイラクにおける2時点の電化率を比較することによって、電力供給の改善状況を検証したいと思います。

2007年と2022年の比較ができれば良いのですが、夜間光データ（月次）は、[DMSP OLS: Nighttime Lights Time Series Version 4](https://developers.google.com/earth-engine/datasets/catalog/NOAA_DMSP-OLS_NIGHTTIME_LIGHTS)（1992年1月～2014年1月）, [VIIRS Nighttime Day/Night Band Composites Version 1](https://developers.google.com/earth-engine/datasets/catalog/NOAA_VIIRS_DNB_MONTHLY_V1_VCMCFG)（2012年4月～2022年6月）のデータに分かれています。

また今回は電化率（電気にアクセスできる人口の割合）を計算するために、[WorldPopの人口分布データ](https://developers.google.com/earth-engine/datasets/catalog/WorldPop_GP_100m_pop)（年次）を使用しますが、イラクの入手できる最新のデータは2020年のものでした。

そのため、今回はVIIRSは2012年8月と2020年8月、人口分布データは2012年と2020年を使用し、2012年と2020年の2時点の電化率を比較することにしました。

# 分析環境
Google Earth Engineを使用しました。


# 手順
1\. 関心領域のインポート（今回はイラクとします）
2\. 夜間光画像のインポート
3\. 人口分布データのインポート
4\. 夜間光の閾値を設定する
5\. 電化率の計算


# 1\. 関心領域のインポート
"country_co"は[ここ](https://en.wikipedia.org/wiki/List_of_FIPS_country_codes)から確認できます。

```js:sample.js
var iraq = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
         .filter(ee.Filter.eq("country_co", "IZ"));

Map.centerObject(iraq, 5);
Map.addLayer(iraq, {}, "iraq");
```

<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/f526d289-04f9-9f9c-eac4-9b388dc93027.png">

# 2\. 夜間光画像のインポート
まずは2012年8月の画像、2020年8月の画像をそれぞれインポートし、表示させます。夜間光は月に1枚しか画像がないため、ee.Imageを使います。clipはインポートした「衛星画像」を上記で設定した境界に合わせて切り取るものです。

```js:sample.js
//2012年8月
var ntl12 = ee.Image("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG/20130801")
              .select("avg_rad")
              .clip(iraq);

//2020年8月
var ntl20 = ee.Image("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG/20200801")
              .select("avg_rad")
              .clip(iraq);

var nighttimeVis = {min: 0.0, max: 100.0};
Map.addLayer(ntl12, nighttimeVis, 'Nighttime light of Iraq in 2012');
Map.addLayer(ntl20, nighttimeVis, 'Nighttime light of Iraq in 2020');
```

2時点の画像を比べると、2020年の方が夜間光の値が大きい（より白く光っている）ように見えます。

<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/8805517a-8f23-3512-5ad3-a3df23c54ce0.png">
<img width="303" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/38c184a3-4340-c24f-a0d0-0c3794028277.png">

左：2012年8月の画像、右：2020年8月の画像
※夜間光の値が大きい（明るい）ほど、白く表示される。

# 3\. 人口分布データのインポート
人口分布データは[WorldPop Global Project Population Data: Estimated Residential Population per 100x100m Grid Square](https://developers.google.com/earth-engine/datasets/catalog/WorldPop_GP_100m_pop) を使用します。解像度は約100mで、2000年1月から2021年1月までのデータが使用可能となっています。

まずは2012年1月～2022年1月までのデータをインポートします。
printによって、Console（コンソール）＞features: List (2241 elements)で、該当期間のデータが全てインポートされていることが確認できます。これを見ると、データは国別に1年に1画像しかありません。また、2021年の画像はなく2020年の画像まで利用可能ということがわかりました。そのため、イラクの2012年と2020年の2つの画像をインポートします。

```js:sample.js
var WorldPop =ee.ImageCollection("WorldPop/GP/100m/pop")
                .filterDate('2012-01-01', '2021-01-01');
 
print("WorldPop:", WorldPop);

//2012年
var pop12 = ee.Image("WorldPop/GP/100m/pop/IRQ_2012")

//2020年
var pop20 = ee.Image("WorldPop/GP/100m/pop/IRQ_2020")

var popVis = {
  bands: ['population'],
  min: 0.0,
  max: 50.0,
  palette: ['24126c', '1fff4f', 'ffff00'] // blue, green, yellow 
};

Map.addLayer(pop12.clip(iraq), popVis, 'Population of Iraq in 2012');
Map.addLayer(pop20.clip(iraq), popVis, 'Population of Iraq in 2020');

```

2時点の画像を比べると、あまり大きな差はないように見えます。

<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/37ac00f8-1d96-3a62-082b-b87299fcb9a2.png">
<img width="296" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/c65109d6-726c-1bf7-a93b-ecad2a9eaa5c.png">

左：2012年の画像、右：2020年の画像
※人口が少ない順に、青（人口が少ない）、緑、黄色（人口が多い）で表示される。

夜間光画像と人口分布画像を比べると、人口が多いところで夜間光の値が大きい（人が集まるところに電力が集中している）ように見えます。

<img width="298" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/38c184a3-4340-c24f-a0d0-0c3794028277.png">
<img width="296" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/c65109d6-726c-1bf7-a93b-ecad2a9eaa5c.png">

左：2020年8月の夜間光画像、右：2020年の人口分布画像

# 4\. 夜間光の閾値を設定する
[論文（Falchetta. et al., 2019）](https://www.nature.com/articles/s41597-019-0122-6)では、ノイズと一時的な光を取り除くために、下限を0.25 nanoWatts/cm2/sr（2014年）、0.35 nanoWatts/cm2/sr（2018年）に設定し、それを超える場合は可視光を検出できる場所（電気にアクセスできる地域）としています。これを踏まえて、今回は2012年8月の夜間光の下限を0.25、2021年8月の夜間光の下限を0.35と設定し、電気にアクセスできる人口を算出します。

まず夜間光の閾値を設定します。閾値以上の地域は1、閾値未満の地域は0という値を割り振りました。
かつ、人がいない地域（湖など）が含まれることを防ぐために、人口が0.1以上である地域（人がいると考えられる地域）に限定しました。

```js:sample.js
//2012
var pop_light_12 = ntl12.gt(0.25).and(pop12.gt(0.1));

//2021
var pop_light_20 = ntl20.gt(0.35).and(pop20.gt(0.1));

var nolightvis = {min:0, max: 1, palette:["black", "white"]};
Map.addLayer(pop_light_12, nolightvis, "Unelectrified area in 2012");
Map.addLayer(pop_light_20, nolightvis, "Unelectrified area in 2021");
```
そして可視光を検出できない場所（電気にアクセスできない地域）を黒色で示したのが下図になります。2012年に比べて、2020年は黒い地域が大幅に減っていることがわかります。

<img width="259" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/a69fbafa-db56-8561-214e-e80e02d1bb34.png">
<img width="254" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2493722/874ec1b7-7943-5f43-5697-c5dc86338674.png">

左：2012年のデータ、右：2020年のデータ
※閾値0.25もしくは0.35未満（電気にアクセスできない地域）の地域は黒、閾値0.25もしくは0.35以上（電気にアクセスできる地域）は白で表示される。

# 5\. 電化率の計算
次に上記で示した白い地域に住む人口（電気にアクセスできる人口）とイラクの総人口を算出します。そしてその電気にアクセスできる人口÷総人口で電化率を計算します。

まず2012年の電化率を計算します。
```js:sample.js
//2012
//Population with access to electricity
var pop_light_12_rev = pop12.updateMask(pop_light_12);

var pop_elect_12 = pop_light_12_rev.reduceRegion({
  reducer: ee.Reducer.sum(),
  scale: 100,
  maxPixels: 1e10
});

print("Population with electricity in 2012:", pop_elect_12);

//Total population
var pop12 = pop12.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: iraq.geometry(),
  scale: 100,
  maxPixels: 1e10
});

print("Total population in 2012:", pop12);
```

printによって、右のコンソールに
Population with electricity in 2012: 26,988,916
Total population in 2012: 29,784,875
が表示されました。

よって、2012年の電化率は
26,988,916÷29,784,875×100 ≒ **90%**
となります。

次に2020年の電化率も計算します。

```js:sample.js
//2020
//Population with access to electricity
var pop_light_20_rev = pop20.updateMask(pop_light_20);

var pop_elect_20 = pop_light_20_rev.reduceRegion({
  reducer: ee.Reducer.sum(),
  scale: 100,
  maxPixels: 1e10
});

print("Population with electricity in 2020:", pop_elect_20);

//Total population
var pop20 = pop20.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: iraq.geometry(),
  scale: 100,
  maxPixels: 1e10
});

print("Total population in 2020:", pop20);
```
printによって、右のコンソールに
Population with electricity in 2020: 37,151,274
Total population in 2020: 38,239,669
が表示されました。

よって、2020年の電化率は
37,151,274÷38,239,669×100 ≒ **97%**
となります。

したがって、2012年の90%に対し2020年は97%となっており、電化率が増加していることがわかりました。

# おわりに
上記の通り、夜間光データを用いて2時点の電化率を計算し、2012年から2020年にかけて電化率が約7%も増加したことがわかりました。ただし、冒頭で示した課題（停電）については、今回の夜間光データでは明らかにできていません。夜間光データで電化率は計算できますが、こうした短期間で変化する電気使用状況については、他のデータや質的な調査によるインタビューの結果から検証する必要があるかと思います。

また、[世界銀行のデータ](https://data.worldbank.org/indicator/EG.ELC.ACCS.ZS?locations=IQ)によると、2012年の電化率は99%、2020年の電化率は100%となっており、若干の乖離があります。国際機関のデータと異なる点が気になりますが、上記の手順により各国レベルだけでなく地域や地方などのより小さい単位での電化率も計算することができるため、ODAプロジェクトのターゲティングや効果検証の際に、参考情報として活用できると考えます。

# 参考文献
Falchetta, G., Pachauri, S., Parkinson, S. et al. A high-resolution gridded dataset to assess electrification in sub-Saharan Africa. Sci Data 6, 110 (2019). https://doi.org/10.1038/s41597-019-0122-6

# コード
上記で示したコードは[こちら](https://code.earthengine.google.com/7d0280a9463cbbb64a2c42ed0193430d)から見れます。