LoginSignup
34
37
@ishimotojuri(Ishimoto Juri)

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

Last updated at Posted at 2023-11-12

はじめに

今回扱うプロジェクト

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

本記事の内容

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

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

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

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

分析方法

分析対象エリア

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

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

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

使用データ

衛星データ

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

地理データ

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

分析方法

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

分析手順

データ収集

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

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

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

データ期間：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:を見ると、どのデータがインポートできたか確認できます。
また、→夜間光データは毎月一枚であることがわかります。
image.png

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

集計・エクスポートのコードはこちら
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」ボタンをクリックすると、画面が表示されます。
image.png

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

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

続けて、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:を見ると、どのデータがインポートできたか確認できます。
また、土地分類データは毎年一枚であることがわかります。
image.png

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

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

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

都市抽出及びエクスポートのコード
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」をクリックします。
image.png

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

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

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

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

  • hotosm_uzb_railways_lines_shp.zip
    image.png

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

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

マップを拡大すると、駅が見つかります。
image.png
image.png
駅を見つけたら、新規ベクタレイヤ（ポイント）を作成します。作成方法はこちらをご参照ください。
image.png

データ作成

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

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

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

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

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

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

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

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

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

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

各メッシュと駅との距離を計算します。
「ベクタ」→「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」とします
image.png
image.png

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

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

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

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

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

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

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

完了後、「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」をダブルクリック→「フィールド属性の編集」→上のアイコン「フィールド計算機」をクリックします。
image.png

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

「レイヤプロパティ」→「単一定義」を「連続値による定義」にします。
image.png
以下の通り設定し、「OK」をクリックします。
値：NL_diff
カラーランプ：Spectral→その上にある「カラーランプを反転」を選択
モード：等量分類
分類数：3
自動で分類後、値をクリックし、以下の通り設定します。
青→上の値：-0.1
黄→上の値：0.1

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

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

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

データのダウンロード

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

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

まとめ

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

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

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

34
37
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
34
37