鉄道は地域経済発展をもたらすの？衛星データで可視化してみた

はじめに

今回扱うプロジェクト

ウズベキスタン国内の鉄道の中でも、カルシ-テルメズ間は、同国から他国を経由せずアフガニスタンへ至る唯一の鉄道路線です。貨物輸送にも重要な役割を果たしており、その鉄道輸送需要はさらなる拡大が予想されています。しかし、カルシ-テルメズ区間には急勾配の山岳地帯が含まれるため、複線化による輸送力増強が困難でした。そこで、同区間において、ディーゼル方式よりも牽引力の高い電気方式への改良を日本が支援しました。


（出所：ODA見える化サイト）

本記事の内容

この支援（事業）によって、対象エリアにおける「経済活性化」が期待されていました。そこで、衛星・地理データを活用して、事業による「経済活性化」を可視化してみたいと思います。

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

最終的なアウトプットは下図です。図では、各指標の変化を色別に表しています。このように図示することによって、どこでどのような変化があったのかがわかりますね。さらに、分析対象駅周辺での夜間光量の変化が、他の地域と比べて大きければ、事業が「経済活性化」に貢献していると推測することができます。
なお、都市については、今回使用したデータによると、2012年から2019年にかけて変化がなかったため、都市と分類された地域だけ黄色で示しています。



（上：夜間光、中央：都市、下：人口）

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

分析方法

分析対象エリア

事業対象エリアはウズベキスタンのカシュカダリア州及びスルハンダリア州（下図の黄色部分）です。
事業で整備した対象線路は同図の赤線部分です。

対象線路に位置する駅（対象駅）において、事業による「経済発展」の影響が顕著に見られると考えるため、分析対象エリアは対象駅とします。
その裏付けとして、先行文献においても、鉄道を利用した地域活性化には、鉄道駅周辺における都市開発が重要だと指摘されています。

また、対象駅は24駅ありますが、分析対象駅として、カルシ駅とテルメズ駅（下図）を選定します。
上記の駅は、事業対象エリアにおいて、その都市人口がそれぞれ最も多い2駅であるため、事業による影響が大きく、事業の貢献度をより検証しやすいと考えるからです。

使用データ

衛星データ

今回の分析には、以下の衛星データを使います。

• 夜間光：VIIRS Nighttime Day/Night Band Composites Version 1

  • データ期間：2012～2023
  • 解像度：約500m

• 土地分類（都市）： MODIS Land Cover Type Yearly Global 500m

  • データ期間：2001～2019
  • 解像度：500m

• 人口：WorldPop: Unconstrained individual countries UN adjusted.

  • データ期間：2000～2020
  • 解像度：1km

地理データ

また、衛星データに加えて、下表に示す地理データも使います。

• 行政区境界：OCHA Uzbekistan - Subnational Administrative Boundaries

  • データ最終更新日：2020年5月

• 線路：HOTOSM Uzbekistan Railways (OpenStreetMap Export) )

  • hotosm_uzb_railways_lines_shp.zip
  • データ最終更新日：2020年7月

• 駅：無料公開データがないため、自作で対応します。

  • QGISにて上記の線路データ及びGoogle Mapを重ね合わせて目視で駅を特定し、ポイントデータを自ら作成します。

分析方法

分析では上記データの解像度に基づき、事業対象エリアにおいて1kmメッシュを作成し、データ抽出・集計を行います。
分析手法としては、各衛星データ（夜間光・都市・人口）の実績値が分析対象駅の一定圏内でどれほど変化したかを、空間的に検証する方法を採ります。

分析手順

データ収集

衛星データの収集
Google Earth Engineを使って、夜間光データを収集します。

まず、分析領域（ポリゴン）を作成します。
一番上の検索欄に「ウズベキスタン」と入れ、マップ画面にウズベキスタンを表示させます。
＋マークをクリックすると、拡大できます。
マップ画面左上に並ぶ■マークをクリックし、ウズベキスタン全体を覆う四角形を描きます。トルクニメスタンやタジキスタンも入ってしまいますが、後で切り取るので問題ありません。

次に、夜間光データをインポートします。
各条件を次のように設定し、下のインポートのコードを実行します。

データ期間：2012-04-01, 2023-06-01（2012年は4月から、2023年は3月まで入手可能）
データ種類：avg_rad
分析領域：geometry（上記のポリゴンで設定した範囲）

インポートのコードはこちら

sample.js

var ntl = ee.ImageCollection(ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG"))
.filterDate('2012-04-01', '2023-03-01')
.select('avg_rad')
.filterBounds(geometry);
            
print(ntl);

「print」は、必要なデータがインポートできたかを確認するためのコードです。
右側の「Console」画面にインポートしたデータの詳細が表示されます。この中のfeatures:を見ると、どのデータがインポートできたか確認できます。
また、→夜間光データは毎月一枚であることがわかります。

最後に、インポートしたデータを集計し、エクスポートします。
夜間光データは各月のため、各年に集計したうえで、エクスポートします。

集計・エクスポートのコードはこちら

sample.js

// 2012
var ntl2012averaged = ntl.filterDate('2012-04-01', '2012-12-31').mean();

Export.image.toDrive({
            image: ntl2012averaged,
            description: 'NighttimeLights_2012_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2013
var ntl2013averaged = ntl.filterDate('2013-01-01', '2013-12-31').mean();

Export.image.toDrive({
            image: ntl2013averaged,
            description: 'NighttimeLights_2013_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2014
var ntl2014averaged = ntl.filterDate('2014-01-01', '2014-12-31').mean();

Export.image.toDrive({
            image: ntl2014averaged,
            description: 'NighttimeLights_2014_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2015
var ntl2015averaged = ntl.filterDate('2015-01-01', '2015-12-31').mean();

Export.image.toDrive({
            image: ntl2015averaged,
            description: 'NighttimeLights_2015_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2016
var ntl2016averaged = ntl.filterDate('2016-01-01', '2016-12-31').mean();

Export.image.toDrive({
            image: ntl2016averaged,
            description: 'NighttimeLights_2016_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2017
var ntl2017averaged = ntl.filterDate('2017-01-01', '2017-12-31').mean();

Export.image.toDrive({
            image: ntl2017averaged,
            description: 'NighttimeLights_2017_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2018
var ntl2018averaged = ntl.filterDate('2018-01-01', '2018-12-31').mean();

Export.image.toDrive({
            image: ntl2018averaged,
            description: 'NighttimeLights_2018_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2019
var ntl2019averaged = ntl.filterDate('2019-01-01', '2019-12-31').mean();

Export.image.toDrive({
            image: ntl2019averaged,
            description: 'NighttimeLights_2019_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2020
var ntl2020averaged = ntl.filterDate('2020-01-01', '2020-12-31').mean();

Export.image.toDrive({
            image: ntl2020averaged,
            description: 'NighttimeLights_2020_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2021
var ntl2021averaged = ntl.filterDate('2021-01-01', '2021-12-31').mean();

Export.image.toDrive({
            image: ntl2021averaged,
            description: 'NighttimeLights_2021_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2022
var ntl2022averaged = ntl.filterDate('2022-01-01', '2022-12-31').mean();

Export.image.toDrive({
            image: ntl2022averaged,
            description: 'NighttimeLights_2022_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2023
var ntl2023averaged = ntl.filterDate('2023-01-01', '2023-3-31').mean();

Export.image.toDrive({
            image: ntl2023averaged,
            description: 'NighttimeLights_2023_Averaged',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

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

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

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

続けて、Google Earth Engineを使って、都市データも収集します。
夜間光データ収集で使ったgeometryを使います。
各条件を次のように設定し、下のインポートのコードを実行します。

データ期間：2012-01-01, 2019-01-01（最新データは2019年1月1日まで入手可能）
データ種類：LC_Type1
分析領域：geometry（上記のポリゴンで設定した範囲）

インポートのコードはこちら

sample.js

var modis = ee.ImageCollection("MODIS/061/MCD12Q1")
            .filterDate("2012-01-01", "2019-12-31")
            .select("LC_Type1")
            .filterBounds(geometry);
            
print(modis);

「print」は、必要なデータがインポートできたかを確認するためのコードです。
右側の「Console」画面にインポートしたデータの詳細が表示されます。この中のfeatures:を見ると、どのデータがインポートできたか確認できます。
また、土地分類データは毎年一枚であることがわかります。

次に、夜間光データと同様に、インポートしたデータをエクスポートします。
そのために、まず都市と分類されたセルだけを抽出し、そのセルを1、それ以外を0とし、データをエクスポートします。

データの詳細ページを開くと、値13は都市ということがわかります。この数字を使って、都市と分類されたセルを抽出します。

以下、上記を実行するためのコードです。

都市抽出及びエクスポートのコード

sample.js

// 2012
var modis_12 = ee.Image("MODIS/061/MCD12Q1/2012_01_01")
var urban_mask_12 = modis_12.eq(13);
var urban_cover_12 = modis_12.updateMask(urban_mask_12);

Export.image.toDrive({
            image: urban_cover_12,
            description: 'Urban_2012',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
});

// 2013
var modis_13 = ee.Image("MODIS/061/MCD12Q1/2013_01_01")
var urban_mask_13 = modis_13.eq(13);
var urban_cover_13 = modis_13.updateMask(urban_mask_13);

Export.image.toDrive({
            image: urban_cover_13,
            description: 'Urban_2013',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
}); 

// 2014
var modis_14 = ee.Image("MODIS/061/MCD12Q1/2014_01_01")
var urban_mask_14 = modis_14.eq(13);
var urban_cover_14 = modis_14.updateMask(urban_mask_14)

Export.image.toDrive({
            image: urban_cover_14,
            description: 'Urban_2014',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
}); 

// 2015
var modis_15 = ee.Image("MODIS/061/MCD12Q1/2015_01_01")
var urban_mask_15 = modis_15.eq(13);
var urban_cover_15 = modis_15.updateMask(urban_mask_15)

Export.image.toDrive({
            image: urban_cover_15,
            description: 'Urban_2015',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
}); 

// 2016
var modis_16 = ee.Image("MODIS/061/MCD12Q1/2016_01_01")
var urban_mask_16 = modis_16.eq(13);
var urban_cover_16 = modis_16.updateMask(urban_mask_16)

Export.image.toDrive({
            image: urban_cover_16,
            description: 'Urban_2016',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
}); 

// 2017
var modis_17 = ee.Image("MODIS/061/MCD12Q1/2017_01_01")
var urban_mask_17 = modis_17.eq(13);
var urban_cover_17 = modis_17.updateMask(urban_mask_17)

Export.image.toDrive({
            image: urban_cover_17,
            description: 'Urban_2017',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
}); 

// 2018
var modis_18 = ee.Image("MODIS/061/MCD12Q1/2018_01_01")
var urban_mask_18 = modis_18.eq(13);
var urban_cover_18 = modis_18.updateMask(urban_mask_18)

Export.image.toDrive({
            image: urban_cover_18,
            description: 'Urban_2018',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
}); 

// 2019
var modis_19 = ee.Image("MODIS/061/MCD12Q1/2019_01_01")
var urban_mask_19 = modis_19.eq(13);
var urban_cover_19 = modis_19.updateMask(urban_mask_19)

Export.image.toDrive({
            image: urban_cover_19,
            description: 'Urban_2019',
            folder: 'Uzbekistan',
            scale: 500,
            region: geometry
}); 

最後に、人口データを収集します。
WorldPop: Unconstrained individual countries UN adjustedのページに行くと、国ごとのデータページが出てくるので、検索欄に「Uzbekistan」と入れて、右の「Data&Resources」をクリックします。

今回は事業計画年（2012）から最新年（2020）までのデータをダウンロードします。
2012年の行にある「Data&Resources」をクリックすると、データのページに飛ぶので、右下の「uzb_ppp_XXXX_1km_Aggregated_UNAdj.tif」をクリックします。
ファイル名は「pop2012」とし、保存します。
これを2013～2020まで繰り返します。

地理データの収集
OCHA Uzbekistan - Subnational Administrative Boundariesのページから、ウズベキスタンの行政区境界データをダウンロードします。
ページの下にある「Data and Resources」から以下のデータをダウンロードします。

• uzb_admbnda_adm2_2018b.zip（地区レベル）→「uzb_district」というファイル名で保存
  

HOTOSM Uzbekistan Railways (OpenStreetMap Export)のページから、ウズベキスタンの線路データをダウンロードします。
ページの下にある「Data and Resources」から以下のデータをダウンロードします。

• hotosm_uzb_railways_lines_shp.zip
  

駅データの作成
駅データについては、無料で公開されているものがないため、QGISで作成します。
QGISを立ち上げ、上記でダウンロードした線路データをQGISのメイン画面にドラッグ＆ドロップして表示させます。
さらに、Google Mapを表示させます。QGISでGoogle Mapを表示させる方法はこちらをご参照ください。

線路データとGoogle Mapを手掛かりに以下の分析対象駅（カルシ駅、テルメズ駅）を探します。
・カルシ駅（カシュカダリア州で人口が最も多い都市）
・テルメズ駅（スルハンダリア州で人口が最も多い都市）

マップを拡大すると、駅が見つかります。


駅を見つけたら、新規ベクタレイヤ（ポイント）を作成します。作成方法はこちらをご参照ください。

データ作成

QGISでプロジェクトCRSを設定します。
メニューの「プロジェクト」→「プロパティ」をクリックします。
上部の「フィルター」欄に「32641」と入力すると、中央の「座標参照系」に
「WGS84/UTM zone 41N」という座標系が表示されるので、それをクリックして「OK」とします。

分析領域を作成します。
上記でダウンロードした「uzb_district.shp」のデータをQGISのメイン画面にドラッグ＆ドロップして表示させます。
ツールバーにある「地物選択」をクリックし、Shiftを押しながら、カシュカダリア州、スルハンダリア州の部分をクリックします。
画面左の「uzb_district」のレイヤを右クリックし、「エクスポート」→「選択地物の保存」をクリックし、「uzb_target_district」というファイル名で保存します。

メッシュデータを作成します。
「ベクタ」→「調査ツール」→「グリッドの作成」で以下の通り設定し、「実行」をクリックします。
グリッドタイプ：長方形（ポリゴン）
グリッドの領域：右側「…」をクリックして「レイヤの領域を使う」でプルダウンから「uzb_target_district」を選択してOK
水平・垂直方向間隔：1000
水平・垂直オーバーレイ：0
Grid CRS：プルダウンから「プロジェクトCRS：EPSG32641…」を選択

「ベクタ」→「空間演算ツール」→「クリップ」で以下の通り設定し、「実行」をクリックします。
入力レイヤ：プルダウンから「グリッドベクタの出力［EPSG：32641］」を選択
クリップレイヤ：プルダウンから「uzb_target_district」を選択

すると、下図のような、分析対象エリアにクリップされたメッシュの図が表示されます。

新たに作成された「出力レイヤ」を右クリック→「エクスポート」→「地物の保存」をクリックし、ファイル名は「mesh.shp」、CRSは「EPSG:32641」で保存します。

メッシュデータと行政区データを統合します。
画面右側の「プロセッシングツールボックス」→「ツールボックス」→「ベクタ一般」→「空間結合（集計つき）」をダブルクリックします。
以下の通り設定し、「実行」をクリックします。
入力レイヤ：mesh
結合するレイヤ：uzb_target_district
ジオメトリの空間関係：「overlaps」と「within」の2つにチェックを入れます
集計する属性：右側「…」をクリックし「ADM0_EN」「ADM1_EN」「ADM2_EN」にチェックを入れて「OK」
計算する集計関数：最大値

新たに作成された「出力レイヤ」を右クリック→「エクスポート」→「地物の保存」をクリックし、ファイル名は「mesh_adm.shp」、CRSは「EPSG:32641」で保存します。

メッシュデータと駅データを統合します。
「station.shp」をQGISのメイン画面にドラッグ＆ドロップします。
「プラグイン」→「プラグインの管理とインストール」で「プラグイン」画面を開きます。

左側メニューで「すべてのプラグイン」が選択された状態で、上部の検索欄に「nnjoin」と入力すると「NNJoin」というプラグインが現れるので、チェックボックスにチェックを入れて「プラグインをインストール」とクリックします。

各メッシュと駅との距離を計算します。
「ベクタ」→「NNJoin」で以下の通り設定し、「OK」をクリックします。
Input vector layer：プルダウンから「mesh_adm」を選び、「Approximate geometries by centroids」にチェックします
Join vector layer：プルダウンから「station」を選び、「join_prefix」を「station」とします
Output layer：新規レイヤ名は「mesh_adm_station」とし、「Neighbour distance field」を「dist」とします

「mesh_adm_station」を右クリックし「属性テーブルを開く」と、「station」との距離（メートル）を示す新たな変数「dist」が作成されていることが分かります。

都市データを追加します
「Urban_2012.tif」をQGISのメイン画面にドラッグ＆ドロップします。
画面右側の「プロセッシングツールボックス」から、「ラスタ解析」→「ゾーン統計量」を選択します。
以下の通り設定し、「実行」をクリックします。
ラスタレイヤ：urban_2012
対象バンド：バンド01:LC_Type1（Gray）
地域ベクタレイヤ：mesh_adm_station
出力カラムの接頭辞：urban12
計算する統計：右側「…」から「合計」のみチェックしてOK
上記を2013～2020も同様に行います。

人口データを追加します。
「pop2012.tif」をQGISのメイン画面にドラッグ＆ドロップします。
画面右側の「プロセッシングツールボックス」から、「ラスタ解析」→「ゾーン統計量」を選択します。

以下の通り設定し、「実行」をクリックします。
ラスタレイヤ：pop2012
対象バンド：バンド1（Gray）
地域ベクタレイヤ：mesh_adm_station
出力カラムの接頭辞：pop12
計算する統計：右側「…」から「合計」のみチェックしてOK
上記を2013～2020も同様に行います。

完了後、「mesh_adm_station」レイヤの属性テーブルを確認すると新たに「pop12sum～pop20sum」という変数が入っていることが確認できます。

夜間光データを追加します。
「NighttimeLights_2012_Averaged.tif」をQGISのメイン画面にドラッグ＆ドロップします。

画面右側の「プロセッシングツールボックス」から、「ラスタ解析」→「ゾーン統計量」を選択します。
以下の通り設定し、「実行」をクリックします。
ラスタレイヤ：NighttimeLights_2012_Averaged
ラスタバンド：バンド1（Gray）
地域ベクタレイヤ：mesh_adm_station
出力カラムの接頭辞：NL12
計算する統計：右側「…」から「平均」のみチェックしてOK
上記を2013～2023も同様に行います。

完了後、「mesh_adm_station」レイヤの属性テーブルを確認すると新たに「NL12mean～NL23mean」という変数が入っていることが確認できます。

都市データを追加します。
「urban_2012.tif」をQGISのメイン画面にドラッグ＆ドロップします。
画面右側の「プロセッシングツールボックス」から、「ラスタ解析」→「ゾーン統計量」を選択します。
以下の通り設定し、「実行」をクリックします。
ラスタレイヤ：urban_2012
ラスタバンド：バンド1（Gray）
地域ベクタレイヤ：mesh_adm_station
出力カラムの接頭辞：urban12
計算する統計：右側「…」から「平均」のみチェックしてOK
上記を2013～2023も同様に行います。

同じレイヤを右クリック→「エクスポート」「地物の保存」でファイル名を「mesh_adm_station」として、上書き保存します。

データ可視化

夜間光データの変化を計算します。
「mesh_adm_station」をダブルクリック→「フィールド属性の編集」→上のアイコン「フィールド計算機」をクリックします。

以下の通り設定し、「OK」をクリックします。
フィールド名：NL_diff
フィールド型：小数点付き数値(real)
フィールド長：10
精度：3
式： "NL23mean" - "NL12mean"
式の隣のボックスの▼フィールドと値の「NL12mean」と「NL23mean」をそれぞれダブルクリックすると、式に反映されます。

「レイヤプロパティ」→「単一定義」を「連続値による定義」にします。

以下の通り設定し、「OK」をクリックします。
値：NL_diff
カラーランプ：Spectral→その上にある「カラーランプを反転」を選択
モード：等量分類
分類数：3
自動で分類後、値をクリックし、以下の通り設定します。
青→上の値：-0.1
黄→上の値：0.1

すると、下図の通り変化が可視化されます。

都市データの変化を計算します。
都市データについても、上記の夜間光と同様の手順で行います。
都市については、今回使用したデータによると、2012年から2019年にかけて変化がありませんでした。

人口データの変化を計算します。
人口データについても、上記の夜間光と同様の手順で行います。

データのダウンロード

最後に、さらにStataで統計分析を行うために、データをダウンロードします。

CSVファイルをエクスポートします。
Stataで更なる分析を行うために、データをCSVファイルの形式でエクスポートします。
「mesh_adm_station」レイヤを右クリック→「エクスポート」→「地物の保存」をクリックします。
以下の通り設定し、「OK」をクリックします。
形式：カンマで区切られた値[CSV]
ファイル名：任意
CSV：EPSG:32641
エンコーディング：UTF-8

まとめ

この記事では、衛星・地理データを活用して、事業による「経済活性化」を可視化してみました。
作成した各指標の変化を図示したところ、分析対象駅の周辺で夜間光と人口が増加しているため、事業が「経済活性化」を促進していると推測できます。
ただ、図からは、詳細な変化の度合いがわかりづらいので、Stata（統計解析ツール）を使って、折れ線グラフ等で詳細を確認したいと思います。
具体的には、分析対象駅から最大50㎞の範囲内で、事業計画時からデータが利用可能な最新年までの各指標の経年変化を分析します。
アウトプットのイメージは下図のとおりです。

この解析方法については、別の記事で紹介します。

ここまで長い記事をご覧いただきありがとうございました！衛星・地理データを活用することで、事業の詳細な影響を確認することができます。
ご関心の分野・地域で、ぜひ試してみてください！質問・コメントもお待ちしております(^^)