見出し画像

夜間光データと人口データを利用した未電化地域の可視化＋電化率の算出の解説 using Google Earth Engine

開発くん

はじめに

先日、夜間光データと人口データを利用してタンザニアの未電化地域を可視化した（＋電化率の算出）記事を投稿したのですが、閲覧数が約100と意外と読まれているようです。

まだ読んでない方は👇

この記事を書いて以降、Google Earth Engine（GEE）に頼ってGISデータ分析しているだけじゃ次のレベルには行けないだろう！！と思い立ち、RでもGISデータ分析が出来るように勉強し始めました。

Rには、rgeeというGEEをR上で動かすパッケージがあるのですが、その習得に多くの時間を割いているという感じです、はい。

まぁそれはいいとして、今回は上の記事で行った分析の解説を軽くしていきます。GEEを使ったGISデータ分析がもっと一般的になればいいなという思いからです。

分析環境

前回と同様、今回使用する分析環境は、我らが最強のGISデータ分析オンラインプラットフォームGEEです。最初は衛星データの処理・分析のためにリリースされたプラットフォームでしたが、後に衛星データ以外のGISデータも追加されました（現在進行形）。

GEEは、Java Scriptと呼ばれるプログラミング言語で動きます。上にも記載したとおり、rgeeというパッケージを用いればRでも動かせます。もちろんのこと、Pythonでも動かすことができます。

本題

と、分析環境を紹介したところで、早速、本題に入っていきます。

分析目的と手順の明確化

データ分析というと、すぐにデータを触りたがる人がいます。けれど、データ分析を行うにあたって重要なことの一つに、自分がどんな目的を持っていて、その目的の達成のためにどういう手順を踏めばいいのかをまず考えることがあると思います（ざくっと方向性だけでも）。

今回の場合、分析目的は、

タンザニアにおける未電化地域の可視化＋電化率の算出

となります。

目的を達成するための手順として以下が考えられます。

  1.  タンザニアの行政区分データをインポート

  2.  夜間光データをインポート

  3.  電化地域と未電化地域の可視化（閾値の設定）

  4.  人口データをインポート

  5.  居住地域と非居住地域の可視化（閾値の設定）

  6.  居住地域を基に未電化地域のみを抽出

  7.  電化率の算出

上記手順に沿って話を進めます。

1. タンザニアの行政区分データをインポート

今回、GEEに格納されている行政区分データを使います。外部データを使ってもいいのですが、なるべくGEE内で完結させた方が再現性が高くなるからです。

コードは至ってシンプルです。

var tz = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
         .filter(ee.Filter.eq("country_co", "TZ"));

これで、タンザニアの行政区分データが適切に抽出できているか確認してみましょう（👈意外と大事）。

Map.centerObject(tz, 6);
Map.addLayer(tz, {}, "Tanzania");
タンザニアの行政区分データの表示

無事に出来ました。次いきましょう、次。

2. 夜間光データをインポート

今回使用する夜間光データは「VIIRS Nighttime Day/Night Band Composites Version 1」（Earth Observation Group, Payne Institute for Public Policy, Colorado School of Minesが提供）となります。

その理由は、2012年4月～のデータが提供されていて、分析対象年の2014年と2020年をカバーしているからです。

では、2014年と2020年の夜間光データをインポートしましょう。コードは以下となります。

var NL2014 = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG")
             .filterDate("2014-01-01", "2015-01-01")
             .select("avg_rad")
             .median()
             .clip(tz);

var NL2020 = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG")
             .filterDate("2020-01-01", "2021-01-01")
             .select("avg_rad")
             .median()
             .clip(tz);

これで、タンザニアにおける2014年と2020年の夜間光データ（中央値）が適切にインポートできたか確認してみましょう（👈意外と大事）。

Map.addLayer(NL2014, {}, "Night Light 2014");
Map.addLayer(NL2020, {}, "Night Light 2020");
左図：2014年の夜間光データ（中央値）
右図：2020年の夜間光データ（中央値）

無事に出来たようです。次いきます。

3. 電化地域と未電化地域の可視化（閾値の設定）

夜間光データを無事にインポートできましたので、これから条件を色々と設定し、電化地域と未電化地域の抽出・可視化に取り組みます。

まず、上で示した夜間光データの図を改めて見ると、夜間光が観測できる場所は白塗り（つまり電化地域）、夜間光が観測できない場所は黒塗り（つまり未電化地域）されていることが分かります。

このデータに対して閾値を設定し、電化地域と未電化地域を切り分けます。

使用している夜間光データは、-1.5～340572.84の値を取ります。値が0を超えると夜間光が観測できているはずなので、その辺りで閾値を設定したいところです。

というわけで、閾値を0.1、0.2、0.3、0.4、0.5の5段階で設定してみて、どの値が閾値として適切かを図で確認してみましょう。

コードは以下となります。関数を使う方法もあるのですが、読者全員が扱えるわけではないと思うので、オーソドックスなコードを記載しました。

var visParams_010 = {min: 0, max: 0.1, palette:["black", "white"]};
var visParams_020 = {min: 0, max: 0.2, palette:["black", "white"]};
var visParams_030 = {min: 0, max: 0.3, palette:["black", "white"]};
var visParams_040 = {min: 0, max: 0.4, palette:["black", "white"]};
var visParams_050 = {min: 0, max: 0.5, palette:["black", "white"]};

Map.addLayer(NL2014, visParams_010, "NL2014_010");
Map.addLayer(NL2014, visParams_020, "NL2014_020");
Map.addLayer(NL2014, visParams_030, "NL2014_030");
Map.addLayer(NL2014, visParams_040, "NL2014_040");
Map.addLayer(NL2014, visParams_050, "NL2014_050");

上記コードを実行すると、地図上で以下の図が表示されるはずです。

上部：左から0.1、0.2、0.3
下部：左から0.4、0.5

上図から、閾値を低く設定しすぎた場合、ノイズのようなものがあり、電化されていない場所まで白くなっていることが分かるかと思います。

なので、丁度良い塩梅の閾値を見定める必要があります（完全に感覚的なものとなります）。

僕個人の印象として、0.4か0.5のどちらかが良さそうを持ちました。ということで、閾値を0.4で設定します。

早速、2014年と2020年の夜間光データに閾値（0.4）を設定します。コードは以下となります（「gt」はgreater thanを意味します）。

var NL2014_gt_040 = NL2014.gt(0.4);
var NL2020_gt_040 = NL2020.gt(0.4);

これで、閾値以下の地域には0、閾値以上の地域には1という値が割り振られました。

上記ができているか、地図で確認していきましょう。0が割り振られた地域は黒色で、1の地域は白色で表示されるようにします。

var visParams_bi = {min: 0, max: 1, palette:["black", "white"]};

Map.addLayer(NL2014_gt_040, visParams_bi, "Night Light 2014 gt 0.4");
Map.addLayer(NL2020_gt_040, visParams_bi, "Night Light 2020 gt 0.4");

以下の地図が表示されるはずです。

左図：2014年の夜間光データ（閾値：0.4）
右図：2020年の夜間光データ（閾値：0.4）

無事に出来ました。良い感じですね。ただ、これで終わりかと思いきや終わりではありません。

というのも、タンザニアの地図を見ていただくといいのですが、北部にはヴィクトリア湖があったり、西部にはタンガニーカ湖があったりと人が住んでいない場所まで黒塗りされています。つまり、未電化地域が過剰にマッピングされてしまっている状態です。そのため、居住地域のみに絞ってマッピングし直す必要があります。

では、どうやって居住地域のみに絞るのか？

答えは、人口データを使う、です。人口データは、人が住んでいる地域にしかマッピングされていないからです。

4. 人口データをインポート

今回使用する人口データは「WorldPop Global Project Population Data: Constrained Estimated Age and Sex Structures of Residential Population per 100x100m Grid Square」（WorldPopが提供）となります。

人口データのインポート方法は少し特殊ですが、以下となります（「216」は人口データ内で割り振られたタンザニアの番号となります）。

var population = ee.Image(ee.ImageCollection("WorldPop/GP/100m/pop_age_sex_cons_unadj")
                 .toList(234)
                 .get(216))
                 .select("population");

無事にインポートできているか地図でも確認した方がいいですが、記事に入れ込むのが面倒になってきたので省きます（すいません、、）。地図に表示するコードは既に載せているので、各自でやっていただければ幸いです。

5. 居住地域と非居住地域の可視化（閾値の設定）

人口データを無事にインポートできたところで、先ほどの夜間光データと同様、居住地域と非居住地域を抽出・可視化するために閾値を設定していきます。

使用している人口データは、0～21171の値を取ります。値が1以上であれば人が住んでいることになりますので、その辺りで閾値を設定します。念のため安全を取って閾値を0.1にします。

コードは以下となります。

var population_gt_010 = population.gt(0.1);

var visParams_pop = {palette:["ff00ff"]};

Map.addLayer(population_gt_010, visParams_pop, "population");

結果は以下のとおりです。

2020年の人口データ（閾値を0.1）

無事に居住地域を可視化できました。これを眺めているだけでも、人は道路沿いに住むんだなぁとか色々と見えてきて、ご飯が食べれそうです（それはないか）。

6. 居住地域を基に未電化地域のみを抽出

以上より、非居住地域を含む未電化地域のデータと居住地域のデータができたので、その2つを重ね合わせると居住地域のみの未電化地域を抽出・可視化することができます。

早速やりましょう。コードは以下となります（表示する図の関係で夜間光の閾値を「less than 0.4」に変えています）。

var NL2020_lt_040 = NL2020.lt(0.4);
var Com = NL2020_lt_040.subtract(population_gt_010);
var pop_unelect = Com.eq(0);
var visParams_unelect = {palette:["black"]};

Map.addLayer(pop_unelect.updateMask(pop_unelect), visParams_unelect, "Population & Unelectrified");
居住地域における未電化地域

というわけで、居住地域における未電化地域の可視化を行いました。上の図がその結果となります。

どうでしょうか。湖などの非居住地域は無事に取り除けて可視化できていることが分かるのではないでしょうか。

これにて、タンザニアおける未電化地域（2020年）の可視化は終わりとなります。

7. 電化率の算出

では最後に、電化率を算出します。コードは以下となります。

var population_rev = population.updateMask(pop_unelect);

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

print("未電化地域の人口", population_unelectrified.get("population"));

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

print("タンザニアの人口",　population_TZ.get("population"));

※人口データ等で閾値を設定すると元々入っていた数値情報が失われるため、1行目で、元の人口データに対して未電化地域の場所情報を掛け合わせることで、未電化地域における人口データを抽出しています。

上記コードを実行すると、GEE画面右上の「Console」上で結果が表示されます。以下の結果が表示されていると思います。

非電化地域の人口：30,946,845人
タンザニアの人口：51,398,080人

以上より、非電化率は約60％（=30,946,845÷51,398,080）と算出することができます。

まとめ

ということで、今回の記事の目的であった、タンザニアおける未電化地域の可視化＋電化率の算出の解説は、無事に達成されました。めでたしめでたしです。

解説をしますと言いつつも、細かいところは省略しています。一から解説となると説明しきれないと思ったからです（現状でさえ約9,000字）。

GEEのコードの説明等は、GEEの公式ガイドに記載されていますので、そちらから適宜学んでいただけますと幸いです。もちろん、その中で分からないことは多々出てくるかと思うのですが、その際は遠慮なくご質問ください。一人でうんうん悩んでいても辛いだけですし。。。

では、これで今回の記事は終わりたいと思います。長文を読んでいただき、ありがとうございました。間違い等ありましたらコメントでご指摘ください（細かい所でこうしておけばよかったなという箇所もあるのですが、仕事でもないので、そのままにしています。。）。

今回のコードのまとめ👇

var tz = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
         .filter(ee.Filter.eq("country_co", "TZ"));

Map.centerObject(tz, 6);
Map.addLayer(tz, {}, "Tanzania");

var NL2014 = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG")
             .filterDate("2014-01-01", "2015-01-01")
             .select("avg_rad")
             .median()
             .clip(tz);

var NL2020 = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG")
             .filterDate("2020-01-01", "2021-01-01")
             .select("avg_rad")
             .median()
             .clip(tz);

Map.addLayer(NL2014, {}, "Night Light 2014");
Map.addLayer(NL2020, {}, "Night Light 2020");

var visParams_010 = {min: 0, max: 0.1, palette:["black", "white"]};
var visParams_020 = {min: 0, max: 0.2, palette:["black", "white"]};
var visParams_030 = {min: 0, max: 0.3, palette:["black", "white"]};
var visParams_040 = {min: 0, max: 0.4, palette:["black", "white"]};
var visParams_050 = {min: 0, max: 0.5, palette:["black", "white"]};

Map.addLayer(NL2014, visParams_010, "NL2014_010");
Map.addLayer(NL2014, visParams_020, "NL2014_020");
Map.addLayer(NL2014, visParams_030, "NL2014_030");
Map.addLayer(NL2014, visParams_040, "NL2014_040");
Map.addLayer(NL2014, visParams_050, "NL2014_050");

var NL2014_gt_040 = NL2014.gt(0.4);
var NL2020_gt_040 = NL2020.gt(0.4);

var visParams_bi = {min: 0, max: 1, palette:["black", "white"]};

Map.addLayer(NL2014_gt_040, visParams_bi, "Night Light 2014 gt 0.4");
Map.addLayer(NL2020_gt_040, visParams_bi, "Night Light 2020 gt 0.4");

var population = ee.Image(ee.ImageCollection("WorldPop/GP/100m/pop_age_sex_cons_unadj")
                 .toList(234)
                 .get(216))
                 .select("population");

var population_gt_010 = population.gt(0.1);
var visParams_pop = {min:0, max: 1, palette:["ff00ff"]};

Map.addLayer(population, visParams_pop, "population");

var NL2020_lt_040 = NL2020.lt(0.4);
var mix = NL2020_lt_040.subtract(population_gt_010);
var pop_unelect = mix.eq(0);
var visParams_unelect = {palette:["black"]};

Map.addLayer(pop_unelect.updateMask(pop_unelect), visParams_unelect, "Population & Unelectrified");

var population_rev = population.updateMask(pop_unelect);

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

print("未電化地域の人口", population_unelectrified.get("population"));

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

print("タンザニアの人口",　population_TZ.get("population"));

いいなと思ったら応援しよう！

開発くん
ブログ・記事をご覧いただき、ありがとうございます！これからも役に立つ情報を発信していきたいと思いますので、応援よろしくお願いします！！