はじめに

このWEBサイトに設置したリンクは，新しいタブ・ウィンドウで開くことを推奨します．
バグの影響で，新しいタブ・ウィンドウでないと開けないリンクが多いためです．
予めご容赦ください．

やること

Rのsfパッケージ及びterraパッケージを用いて，空間分析を行う方法を紹介します．同じ筆者による『事例で学ぶ経済・政策分析のためのGIS入門』のためのR入門：空間解析・可視化編や『事例で学ぶ経済・政策分析のためのGIS入門』のためのR入門：空間解析・可視化編【補足】，3-5_高速道路開業と雇用水準の変化では，どちらかと言うと学術研究に用いられそうなテクニックを紹介しています．一方このページは，実務上直面しがちな課題への取り組み方について解説します．
　各セクションに共通する流れは，「RとGISできる系の人」という認識を本当はそんなにできないのに社内で持たれている私が，上司から「これ調べられるでしょ？」とメールで無茶振りを受け，それに答えるためのソリューションをRで実装するというものです．事例1はベクタデータ，事例2・3はラスタデータを用いた分析を行います．
　このページの内容は，Rに関する基礎知識を前提として記述されています．もし初めてRに触れる場合，WEBページとしては，卒業論文のためのR入門 (T. Mori)や，『事例で学ぶ経済・政策分析のためのGIS入門』のためのR入門：基礎編などで事前学習を行うことをお勧めします．書籍であれば，『RユーザのためのRStudio［実践］入門』の3・4章がお勧めです．
　なお，このページで用いられるデータの前処理・原典については，その無茶振り，（Rで）GISが解決します：付録を参照してください．加えて，このWEBページの元になっているRmd（R_de_GIS_v10.Rmd）と演習データは，こちらのDropboxリンクに設置してあります．

GISの基礎

GISに関する基礎知識を身につけたい場合，以下の書籍が参考になります．『地図リテラシー入門』は，GISを扱う際の諸注意を分かりやすく解説しています．『事例で学ぶ』の主眼は経済・政策分析にありますが，同書の1章では，実務での活用の際にも共通して必要とされるGISの基礎知識を端的にまとめています．『GISのビジネス活用』は，GISの基礎知識と共に，実務における様々な場面でのGISの活用事例を紹介しています．

　GISを用いて何かしらの実践的課題の解決を行うという，このページのコンセプトは，以下の書籍にインスパイアされたものです．日常生活やビジネス，まちづくりで直面しがちな問題に対して，(Q)GISを用いたソリューションを与える文献です．

その他のGIS関連のチュートリアルのまとめサイトとして，GIS: Useful Links (Mizuki Kawabata Lab)を参照することもお勧めします．

更なる学習のために

sfパッケージとdplyrパッケージを組み合わせたGISに関して，さらに詳細・高度なジオプロセシングを解説しているWEBサイトをいくつか紹介します．

日本語のサイトとしては，

英語のサイトとしては，以下がおすすめです．

アップデート

  • 国土数値情報が提供しているデータのURL変更に対応しました．また，ラスタ解析のセクションを追加しました（2020年6月23日）．
  • 筆者WEBページの移動に伴う各種リンクのURL変更に対応しました．また，外部チュートリアルサイトのまとめを更新しました（2020年8月7日）．
  • 一部データ・sfパッケージの仕様変更に対応しました．また，ネットワーク解析のセクションを追加しました（2020年12月7日）．
  • 『事例で学ぶ経済・政策分析のためのGIS入門』の刊行に伴い，コンテンツを大幅に変更しました（2022年12月18日）．
  • OSMからのデータ取得方法に仕様変更があったため，付録のコードを改訂しました．それに伴い，OSMから取得された店舗ポイントデータも最新のものに差し替えられています（2025年4月25日）．

事前準備

使うパッケージを読み込みます．インストールが済んでいない場合は，予めインストールを済ませておいてください．インストールの過程で「Do you want to install from sources the packages which need compilation? (Yes/no/cancel)」というメッセージが表示された場合には，とりあえず「no」を選択することをお勧めします．

  • tidyverse：各種データ操作のためのパッケージ
  • sf：空間データ操作のためのパッケージ
  • terra：空間データ（ラスタ）操作の為のパッケージ
  • openxlsx：Rでxlsx形式のファイルを読み込み・書き出しするためのパッケージ
  • mapview：インタラクティブな空間データの可視化を行うためのパッケージ
  • RColorBrewer：様々なプロットで使えるカラーパレットを収録したパッケージ
  • osrm：OpenStreetMap上の道路データを用いて経路探索を行うパッケージ
  • movecost：コスト最小経路解析を行うためのパッケージ
library(tidyverse)
library(sf)
library(terra)
library(openxlsx)
library(mapview)
#mapviewパッケージ使用時のおまじない
mapview::mapviewOptions(fgb=FALSE)
library(RColorBrewer)
library(osrm)
library(movecost)

動作環境によっては，sfパッケージの読み込み時にエラーが出る場合があります．その原因は，sfパッケージが依存しているRcppのバージョンが古いことが多いようです．もしエラーが出る場合には，以下のコードを実行して見てください（筆者のPCではエラーが出なかったので，コメント化してあります）．

# install.packages("Rcpp")
# library(Rcpp)

加えて，osrmのいくつかの関数では，バージョンによってはバグによるエラーが出る状態です．それらは最新バージョンでは対応済みとのことなので，もしエラーが出る場合には，以下のコードを実行して見てください（筆者のPCではエラーが出なかったので，コメント化してあります）．

# remotes::install_github("riatelab/osrm")
# library(osrm)

もし，このRmdをクリックしてRStudioを起動する場合は，Rmdが設置された場所が作業ディレクトリ（ざっくり言うと，データの読み込み・書き出し先）として自動的に指定されます．一方，手動で作業ディレクトリを指定する場合には，setwd関数を用います．
　例えば「C:/ktakano/econome」というパスに設定したい場合，以下のように書きます．架空のディレクトリにつき，そのままコードとして実行してもエラーが出るので，ここではコメント化しています．以下では，全てのデータがRmdと同じ階層のフォルダに設置されている前提で，説明を行います．

# setwd("C:/ktakano/econome")

事例1：移動販売事業

今回の無茶振り

大手リテール企業に勤める私に，上司からこのようなメールが届きました．

開発後数十年が経過したニュータウンでは，買い物難民の問題が深刻化しつつあります．
そこで，我が社の新事業として，買い物難民が発生している（可能性がある）地域を対象に，移動販売を展開したいと考えています．
つきましては，部署の皆さんに以下の2つ事項を依頼します．

・買い物難民が発生していそうなニュータウンとその特徴の調査・特定
・それらニュータウンを巡回する最短経路の調査

　買い物難民問題とは，公共交通や自家用車等の移動手段が使えなくなることで，生活必需品の買い物が困難になった地区が，都市内に点在する問題です．こうした問題を打開するため，あるいはそれを新たなビジネスチャンスとして捉えて，様々な大手スーパーが，近年移動販売事業を展開しています．
　事業展開の対象地域は福岡県福岡市で，地域単位は500mメッシュ（500m四方の領域）です．買い物難民地域は差し当たり「市内の平均よりも店舗へのアクセス性が低く，高齢世帯割合が高い地域」と定義します．
　アクセス性の評価方法は様々ありますが，今回は500mメッシュの重心点から1km圏内を各メッシュの商圏とし，その中にある店舗の数で計測します．例えば，Figure 1の例では，アクセス性の値は4です．

Figure 1
Figure 1：店舗アクセスの考え方

データの読み込み

上司からのリクエストに応えるためには，地域の世帯属性や店舗アクセス性，ニュータウンの位置に関する情報が必要そうです．それぞれに対応するデータとして，以下に挙げたものを用意しました．加えて，買い物難民の発生には交通アクセス性が大いに関係しているので，鉄道駅の可用性に着目して，地域属性のひとつとして計測します．

ベクタデータ

今回用意したベクタデータは全てGeoJSON形式です．GeoJSONの読み込みにはsfのread_sf関数を用います．今回は，筆者が用意したDropboxリンクから直接ダウンロードします．

#小地域（福岡市）
chome <- sf::read_sf(dsn="https://www.dropbox.com/s/5614iigqh5nmhtb/chome_fukuoka_c.geojson?dl=1")
#500mメッシュ（福岡市）
mesh <- sf::read_sf(dsn="https://www.dropbox.com/s/lnr54m91keqbyp4/mesh_fukuoka_0500_c.geojson?dl=1")
#店舗（福岡市周辺）
shop <- sf::read_sf(dsn="https://www.dropbox.com/s/a3bcakej2xez65s/shop_fuk_p.geojson?dl=1")
#ニュータウン（福岡県）
nt <- sf::read_sf(dsn="https://www.dropbox.com/s/wmq2psw2k87ef1m/newTown_fukuoka.geojson?dl=1")
#鉄道駅（全国）
sta <- sf::read_sf(dsn="https://www.dropbox.com/s/l6dt8rckuwu2e8a/rail_jp.geojson?dl=1")

読み込んだデータのうち，例えば小地域ポリゴンchomeの先頭行を表示してみます．表示される各情報の意味は以下の通りです．

  • Geometry type：ベクタデータの種類（例：ポリゴン，ポイント，ライン）．
  • Dimension：データの次元．3DのGISデータを読み込んだ場合，XYZとなる．
  • Bounding box：データ内の全地物を覆う最小の長方形の頂点座標．
  • Geodetic CRS：データに定義された座標系．

地物の幾何的属性をWKT（Well-known text）形式で格納した変数がgeometryで，それ以外の変数は各地物の属性を表します．

head(chome)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -484269.7 ymin: 3769765 xmax: -482761 ymax: 3771980
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 5
##   KEN_NAME GST_NAME CSS_NAME MOJI                                       geometry
##   <chr>    <chr>    <chr>    <chr>                                 <POLYGON [m]>
## 1 福岡県   福岡市   東区     箱崎１丁目 ((-483050.8 3769971, -483054.8 3769957,…
## 2 福岡県   福岡市   東区     箱崎２丁目 ((-483491.9 3770205, -483531.1 3770159,…
## 3 福岡県   福岡市   東区     箱崎３丁目 ((-483207.9 3770558, -483208.8 3770559,…
## 4 福岡県   福岡市   東区     箱崎４丁目 ((-483840 3770457, -483848 3770446, -48…
## 5 福岡県   福岡市   東区     箱崎５丁目 ((-483432.5 3771149, -483438.5 3771136,…
## 6 福岡県   福岡市   東区     箱崎６丁目 ((-483384.7 3770987, -483485.6 3771032,…

空間分析の際には，shapefile形式のデータもよく用いられます．例えば，Rmdと同じフォルダ（あるいは指定した作業ディレクトリ）にあるニュータウンデータを読み込みたいときには，同じくread_sf関数を用いて読み込みます．GeoJSONのデータをローカルから読み込む時も同様です．
　ただし，先ほどに加え，optionsという引数を指定しています．「ENCODING=CP932」は，「データの文字コードをCP932（Windowsで作られたデータに多い）にする」という意味です．
　文字コードを指定しない場合，Macで読み込もうとすると文字化けが起こる場合があります．一方，文字コードがUTF-8のデータをWindowsで読み込もうとすると，文字化けが起こる場合があります．
故に，日本語データの読み込みで文字化けが起きる場合には，CP932かUTF-8を試してみると良いと思います．

#ニュータウン（福岡県）
nt <- sf::read_sf(dsn="newTown_fukuoka.shp",
                  options="ENCODING=CP932")

　GeoJSON形式は1つのファイルで空間データとして完結しているため，shapefileに比してハンドリングを行いやすい等の利点があります．一方，shapefile形式が空間データのデファクト・スタンダードである現状は当分変わらなさそうなので，両方の扱い方を知っておくことが望ましいです．

テキストデータ

　今回用いるテキスト形式のデータは，csv形式及びxlsx形式です．csvはread.csv関数，xlsxはopenxlsxのread.xlsx関数で読み込みます．csvの読み込みの際にも，場合によっては引数fileEncodingで，文字コードを指定する必要があります．

#国勢調査メッシュ集計（福岡市）
census <- openxlsx::read.xlsx(xlsxFile="https://www.dropbox.com/s/408xzo8lo6x1idl/census_2015.xlsx?dl=1")
#店舗（福岡市周辺）
shop_tab <- read.csv(file="https://www.dropbox.com/s/jkfxncnjvz5re6g/fuk_p_lonlat.csv?dl=1",
                     fileEncoding="UTF-8")

投影法（座標系）に関する設定

今回の演習データは，全て私の方で座標系を揃えていますが，一般に利用可能な空間データは各々異なる座標系の下で作成されています．本来は，それらをひとつに統一する作業（投影変換と呼ばれます）が必要です．参考までに，現状は投影座標系（WGS84/UTM54N）のニュータウンポイントを，sfのst_transform関数を用いて，地理座標系（WGS84）に投影変換する方法を示します．

nt_lonlat <- nt %>%
  sf::st_transform(crs=sf::st_crs(4326))
head(nt_lonlat)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 130.6879 ymin: 33.8085 xmax: 130.9045 ymax: 33.90049
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 5
##   pref   city     ntname                   year_c            geometry
##   <chr>  <chr>    <chr>                     <int>         <POINT [°]>
## 1 福岡県 北九州市 ユーフォリアヒルズ舞ヶ丘   1985 (130.9045 33.82573)
## 2 福岡県 北九州市 永犬丸・則松               1995 (130.7295 33.84641)
## 3 福岡県 北九州市 下上津役永犬丸             1979 (130.7391 33.83491)
## 4 福岡県 北九州市 高須                       1972 (130.6915 33.89044)
## 5 福岡県 北九州市 若松西部                   1986 (130.6879 33.90049)
## 6 福岡県 北九州市 千代ニュータウン           1976  (130.7429 33.8085)

「crs=sf::st_crs(4326)」が，座標系を指定している箇所になります．4326は，EPSGコードと呼ばれる形式で付与されたWGS84のIDです．任意の地域について，適当な座標系を調べる為のサービスとして，epsg.ioがあります．加えて，日本の空間データでよく用いられる地理・投影座標系についてのEPSGコードは，EPSGコード一覧表/日本でよく利用される空間座標系（座標参照系）で分かりやすくまとめられています．

とりあえず可視化

読み込んだ空間データをとりあえず可視化したい時には，mapviewのmapview関数を用います．例として，店舗ポイントshopを可視化します．背景地図は，インターネットに接続した状態でのみ利用可能です．

mapview::mapview(x=shop)

複数の空間データを重ねる際や，塗りつぶしの色を変えたい場合，引数col.regionを以下のように指定します．例として，小地域ポリゴン（灰）の上に，店舗ポイント（赤）を重ねて可視化します．

mapview::mapview(x=chome,col.region="gray") +
  mapview::mapview(x=shop,col.region="red")

テキストデータの結合・可視化

500mメッシュポリゴンに，テキスト形式の国勢調査データを結合し，変数を可視化します．

国勢調査データの加工

国勢調査データ上で，いくつかの変数を作成します．まず，以後の分析で扱いやすいように，世帯総数をhhという変数に複製します．加えて，高齢単身の一般世帯数と高齢夫婦のみの一般世帯数の和（高齢世帯総数）を世帯総数で割ることで，高齢世帯eld_hh_rateを作成します．
　その後，メッシュコードKEY_CODEと複製・作成した変数のみを残した上で，KEY_CODEを文字列型に変換します．メッシュポリゴン側のKEY_CODEが文字列型のため，そちらに変数型を合わせます．

census_sum <- census %>%
  #変数の作成
  dplyr::mutate(hh=世帯総数,
                #高齢世帯割合
                eld_hh_rate=(高齢単身の一般世帯数+高齢夫婦のみの一般世帯数)/世帯総数) %>%
  #メッシュコードと作成した変数だけ残す
  dplyr::select(KEY_CODE,hh,eld_hh_rate) %>%
  #メッシュコードを文字列型に変換
  dplyr::mutate(KEY_CODE=as.character(KEY_CODE))

500mメッシュに国勢調査を結合

メッシュポリゴンに，上で作成した国勢調査の変数を結合します．

mesh_census <- mesh %>%
  #国勢調査変数を，メッシュコードをキーに結合
  dplyr::left_join(y=census_sum,by="KEY_CODE") %>%
  #変数が結合できたメッシュだけ残す
  dplyr::filter(!is.na(hh))

各変数に基づいて可視化

国勢調査変数に基づいて，メッシュを色分けして可視化します．特定の変数について色分けを行う際には，mapview関数内の引数zcolで，変数名を指定します．例として，世帯総数hhで色分けします．

#世帯総数（hh）の可視化
mapview::mapview(x=mesh_census,zcol="hh")

RColorBrewerのbrewer.pal関数を用いて，塗りつぶしの色分けに使うパレットを変えることができます．以下の例は，青から赤に変わるパレットを使って世帯総数hhで色分けを行っています．また，メッシュの枠線の色は引数colorで指定できます．今回は透明で指定しました．

#世帯総数（hh）の可視化
mapview::mapview(x=mesh_census,
                 zcol="hh",
                 #青から赤に変わるパレットを使って色分け
                 col.regions=rev(RColorBrewer::brewer.pal(n=8,name="Spectral")),
                 #メッシュの枠線を透明にする
                 color="transparent")

同じ手順で，高齢世帯割合eld_hh_rateで色分けした地図も作成してみましょう．

#高齢世帯割合（eld_hh_rate）の可視化
mapview::mapview(x=mesh_census,
                 zcol="eld_hh_rate",
                 col.regions=rev(RColorBrewer::brewer.pal(n=8,name="Spectral")),
                 color="transparent")

メッシュ毎の商圏の定義

各メッシュの商圏を，メッシュの重心から1km圏内の領域と定義します．このような領域は，バッファという機能を用いて生成できます．

メッシュの重心ポイントを生成

重心の生成には，st_centroid関数を用います．

mesh_census_cent <- sf::st_centroid(x=mesh_census)
head(mesh_census_cent)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -493723.4 ymin: 3753364 xmax: -490456.7 ymax: 3756053
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 4
##   KEY_CODE             geometry    hh eld_hh_rate
##   <chr>             <POINT [m]> <dbl>       <dbl>
## 1 503012591 (-490456.7 3753364)     5       0.4  
## 2 503012674 (-492066.2 3754941)    19       0.368
## 3 503012681 (-491529.7 3754415)    51       0.412
## 4 503012682 (-490945.1 3754355)     5       0.4  
## 5 503012763 (-493723.4 3756053)    11       0.636
## 6 503012771 (-492602.5 3755467)     5       0.4

重心ポイントから1kmバッファを生成

バッファの生成には，st_buffer関数を用います．引数distでバッファの距離を指定します．現状の座標系だと距離はm単位なので，1000mとします．

mesh_census_cent_bf <- sf::st_buffer(x=mesh_census_cent,
                                     dist=1000)
head(mesh_census_cent_bf)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -494723.4 ymin: 3752364 xmax: -489456.7 ymax: 3757053
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 4
##   KEY_CODE                                            geometry    hh eld_hh_rate
##   <chr>                                          <POLYGON [m]> <dbl>       <dbl>
## 1 503012591 ((-489456.7 3753364, -489458.1 3753312, -489462.2…     5       0.4  
## 2 503012674 ((-491066.2 3754941, -491067.5 3754889, -491071.6…    19       0.368
## 3 503012681 ((-490529.7 3754415, -490531.1 3754363, -490535.2…    51       0.412
## 4 503012682 ((-489945.1 3754355, -489946.5 3754303, -489950.6…     5       0.4  
## 5 503012763 ((-492723.4 3756053, -492724.8 3756001, -492728.9…    11       0.636
## 6 503012771 ((-491602.5 3755467, -491603.9 3755415, -491608 3…     5       0.4

バッファの可視化

全て可視化すると時間がかかる（最悪Rがクラッシュする）ので，重心・バッファのうち最初の6個を可視化します．

#重心の可視化
mapview(x=head(mesh_census_cent),
        col.region="red") +
  #1kmバッファの可視化
  mapview(x=head(mesh_census_cent_bf))

店舗データの分析

日用品が買える店に絞る

店舗ポイントを，日用品が買える店舗のみに絞ります．

shop_d <- shop %>%
  dplyr::filter(shop%in%c("mall","convenience","supermarket","department_store"))

重心バッファと日用品店ポイントの共通部分を取る

メッシュ重心バッファと日用品店ポイントの空間的な共通部分を取った上で，それぞれの属性を付与します．こうした操作は，sfのst_intersection関数で可能です．これにより，各重心バッファに含まれる駅が特定できます．

bf_shop <- sf::st_intersection(x=mesh_census_cent_bf,
                               y=shop_d)

メッシュ毎に店舗数を集計

重心バッファ毎に，それに含まれる店舗数を集計します．st_drop_geometry関数を用いてgeometry属性（空間データとしての情報）を削除した上で，通常の表形式データと同じ方法で集計作業が可能です．
まず，dplyrのgroup_by関数を用いてメッシュコードKEY_CODE（すなわち重心バッファのID）毎にレコードをグループ化し，summarise関数を用いて，レコード数をカウントします．

bf_shop_sum <- bf_shop %>%
  #geometry属性を削除
  sf::st_drop_geometry() %>%
  #メッシュコード毎にレコードをグループ化
  dplyr::group_by(KEY_CODE) %>%
  #レコード数をカウント
  dplyr::summarise(shop_d=n()) %>%
  #グループ化を解除
  dplyr::ungroup()

国勢調査メッシュに集計結果を結合

国勢調査メッシュに，上で集計した店舗数を結合します．

mesh_census_dshop <- mesh_census %>%
  dplyr::left_join(y=bf_shop_sum,by="KEY_CODE") %>%
  #店舗数がNA（結合できなかった）場合は0に置換
  dplyr::mutate(shop_d=ifelse(!is.na(shop_d),shop_d,0))

店舗数を可視化

結合された店舗数を，メッシュ毎に可視化します．博多・天神地区のような，福岡市の中心市街地では店舗数が多く（店舗アクセス性が高く），そこから離れるほど少なくなることが見て取れます．

mapview(x=mesh_census_dshop,
        zcol="shop_d",
        col.regions=rev(RColorBrewer::brewer.pal(n=8,name="Spectral")),
        color="transparent")

買い物難民化リスクの評価

高齢世帯率・店の多さの2変数から，各メッシュの買い物難民化リスクを算定します．

分類変数の作成

高齢世帯率と店の多さの2変数に基づき，メッシュを分類します．分類変数の名前をaccess_cとします．その定義は以下の通りです．

  • eld_ma：高齢世帯率が平均より高い＆店舗数が平均より多い
  • eld_la：高齢世帯率が平均より高い＆店舗数が平均より少ない
  • young_ma：高齢世帯率が平均より低い＆店舗数が平均より多い
  • young_la：高齢世帯率が平均より低い＆店舗数が平均より少ない

eld_laに該当するメッシュが，潜在的に買い物難民が発生しやすい場所になると考えられます．

mesh_census_dshop_ref <- mesh_census_dshop %>%
  #変数access_cの追加
  dplyr::mutate(access_c=case_when(
    (eld_hh_rate>mean(eld_hh_rate))&(shop_d>mean(shop_d)) ~ "eld_ma",
    (eld_hh_rate>mean(eld_hh_rate))&(shop_d<mean(shop_d)) ~ "eld_la",
    (eld_hh_rate<mean(eld_hh_rate))&(shop_d>mean(shop_d)) ~ "young_ma",
    (eld_hh_rate<mean(eld_hh_rate))&(shop_d<mean(shop_d)) ~ "young_la",
  ))

分類変数が正しく作成できたか，高齢世帯率eld_hh_rateと店舗数shop_dから散布図を書くことにより，確認します．散布図の作成には，ggplot2のgeom_point関数を用います．

#x軸にeld_hh_rate，y軸にshop_dを取り，点をaccess_cで色分け
plot_mesh_census_dshop_ref <- ggplot2::ggplot(data=mesh_census_dshop_ref,
                                              aes(x=eld_hh_rate,y=shop_d,color=access_c)) +
  #散布図を作成
  ggplot2::geom_point() +
  #eld_hh_rateの平均値をx軸方向に点線で描画
  ggplot2::geom_vline(xintercept=mean(mesh_census_dshop_ref$eld_hh_rate),lty=2) +
  #shop_dの平均値をy軸方向に点線で描画
  ggplot2::geom_hline(yintercept=mean(mesh_census_dshop_ref$shop_d),lty=2)
#可視化
plot_mesh_census_dshop_ref

分類変数に応じた可視化

分類変数access_cに基づいて，メッシュを可視化します．

mapview::mapview(x=mesh_census_dshop_ref,
                 zcol="access_c",
                 color="transparent")

買い物難民化リスクが高いニュータウンの特定

ニュータウンポイントにメッシュ（分類変数結合済）を結合

ニュータウンポイントに，メッシュ上で計算された分類変数を空間結合します．

nt_dshop_ref <- sf::st_intersection(x=nt,
                                    y=mesh_census_dshop_ref)

高齢世帯が多く店舗が少ないニュータウンを抽出

高齢世帯が多く店舗が少ないニュータウンポイント（access_cが「eld_la」に一致）を抽出します．

nt_dshop_ref_eld_la <- nt_dshop_ref %>%
  dplyr::filter(access_c=="eld_la")

メッシュとニュータウンを同時に可視化

メッシュとニュータウンポイントをそれぞれ分類変数で色分けし，同時に可視化します．

mapview::mapview(x=mesh_census_dshop_ref,
                 zcol="access_c",
                 color="transparent",
                 #凡例の名前を手動で設定
                 layer.name="Age-commercial access") +
  mapview::mapview(x=nt_dshop_ref_eld_la,
                   col.region="black")

ニュータウンの交通アクセス性

ここまでの分析で，世帯属性と店舗アクセス性の2つの変数から，買い物難民リスクが高いニュータウンを特定しました．一方，これらニュータウンは，実際に公共交通不便地に位置しているでしょうか？鉄道駅へのアクセス性を計測する事で，このことを確かめましょう．
本来であればバス等も考慮した分析が望ましいのですが，分析が複雑になりそうなので，ここでは「鉄道駅を中心としてバス網が広がる」という，日本の公共交通にありがちな構造を想定した上で分析を行います．

メッシュ重心バッファと駅ポイントの共通部分を取る

メッシュ重心バッファと駅ポイントの共通部分を取ります．これにより，各重心バッファに含まれる駅が特定できます．

bf_sta <- sf::st_intersection(x=mesh_census_cent_bf,
                              y=sta)

メッシュコード（KEY_CODE）毎に共通部分の個数を集計

店舗ポイントの場合と同じ要領で，重心バッファ毎に含まれる駅の個数を集計します．

bf_sta_sum <- bf_sta %>%
  #geometry属性を削除
  sf::st_drop_geometry() %>%
  #データをメッシュコードでグループ化
  dplyr::group_by(KEY_CODE) %>%
  #グループ毎にレコード数を集計
  dplyr::summarise(nsta=n()) %>%
  dplyr::ungroup()

集計をメッシュ（分類変数結合済み）に結合

集計された到達可能駅数を，メッシュに結合します．

mesh_census_dshop_ref_t_acc <- mesh_census_dshop_ref %>%
  #メッシュコードをキーに，上の集計結果（到達可能駅数）を結合
  dplyr::left_join(y=bf_sta_sum,by="KEY_CODE") %>%
  #駅数がNAの場合，0に置き換える
  dplyr::mutate(nsta=ifelse(!is.na(nsta),nsta,0))

利用可能駅の個数の可視化

メッシュに結合された，利用可能駅数を可視化します．色分けの設定は，「各変数に基づいて可視化」セクションで使ったものを流用します．

mapview::mapview(x=mesh_census_dshop_ref_t_acc,
                 zcol="nsta",
                 col.regions=rev(RColorBrewer::brewer.pal(n=8,name="Spectral")),
                 color="transparent")

ニュータウンと一緒に可視化

先ほど特定された，買い物難民リスクが高い（高齢世帯が多く店舗が少ない）ニュータウンポイントと一緒に，利用可能駅数を可視化します．

mapview::mapview(x=mesh_census_dshop_ref_t_acc,
                 zcol="nsta",
                 col.regions=rev(RColorBrewer::brewer.pal(n=8,name="Spectral")),
                 color="transparent") +
    mapview::mapview(x=nt_dshop_ref_eld_la,
                     col.region="black")

巡回経路の探索

osrmのosrmTrip関数を用いて，買い物難民リスクが高いニュータウンポイントを巡回する経路の探索を行います．経路探索の際に用いられるのは，巡回セールスマン問題（traveling salesman problem; TSP）と呼ばれる数理最適化のアルゴリズムです．

#最短経路の探索
tsp_eld_la <- osrm::osrmTrip(loc=nt_dshop_ref_eld_la,
                             returnclass="sf")
#返り値の表示
tsp_eld_la
## [[1]]
## [[1]]$trip
## Simple feature collection with 8 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -500676.6 ymin: 3760910 xmax: -478358.7 ymax: 3779731
## Projected CRS: WGS 84 / UTM zone 54N
##   start end  duration distance                       geometry
## 1     1   3 16.045000  18.4378 LINESTRING (-500582.4 37677...
## 2     3   4 20.216667  16.9941 LINESTRING (-486481.7 37613...
## 3     4   7 10.686667   7.4898 LINESTRING (-478358.7 37732...
## 4     7   8  4.101667   3.1036 LINESTRING (-480418.8 37781...
## 5     8   5  8.121667   6.9845 LINESTRING (-481811.7 37796...
## 6     5   6  1.008333   0.3633 LINESTRING (-482291.5 37746...
## 7     6   2 17.008333  20.7410 LINESTRING (-482258.9 37744...
## 8     2   1  9.425000   8.9630 LINESTRING (-495557.1 37645...
## 
## [[1]]$summary
## [[1]]$summary$duration
## [1] 86.61333
## 
## [[1]]$summary$distance
## [1] 83.0771

経路探索の結果はlist形式で返され，1つ目の要素には経路のラインデータ，2つ目の要素には総所要時間[分]，3つ目の要素には総移動距離[km]が入っています．ここでは，1つ目の要素を取り出した上で，ニュータウンポイントと共に地図上にプロットします．

#ニュータウンをプロット
mapview::mapview(x=nt_dshop_ref_eld_la) +
  #探索された経路を黒でプロット
  mapview::mapview(x=tsp_eld_la[[1]]$trip,color="black")

ここまでで，ようやく上司のリクエストに答えられそうな最低限の資料が作れた感じがします．

実践問題

トラックが2台使える（あるいは2回に分けて巡回する日程が組める）ことになったので，福岡市の東地区（中央区，博多区，東区）と西地区（それ以外の区），それぞれ分けて巡回を行うプランを立てることにします．この時，各地区について別個に経路探索を行い，その結果を可視化してください．また，総走行距離は一度に巡回する場合に比べてどう変わったでしょうか？

【分析の準備】行政区ポリゴンの作成

小地域ポリゴンを融合して行政区ポリゴンを作成した上で，福岡市を2つの地区に分けることを考えます．ポリゴンの融合は，dplyrのgroup_by関数とsummarise関数を用いて実行できます．
　小地域ポリゴン内の区名変数はCSS_NAMEなので，まずはCSS_NAMEで小地域をグループ化し，summarise関数でCSS_NAMEの値毎に小地域を融合します．最後に，融合されたポリゴンに対して地区区分変数を追加します．

ward <- chome %>%
  #区名変数でポリゴンをグループ化
  dplyr::group_by(CSS_NAME) %>%
  #ポリゴンをまとめる
  dplyr::summarise() %>%
  dplyr::ungroup() %>%
  #まとめたポリゴンに対し，地区区分変数を追加
  dplyr::mutate(region=case_when(
    CSS_NAME%in%c("中央区","博多区","東区") ~ "東地区",
    TRUE ~ "西地区"
  ))

作成したポリゴンを，地区区分毎に色分けして可視化すると以下のようになります．

mapview::mapview(x=ward,
                 #色分けの基準は地区区分
                 zcol="region")

各ニュータウン（高齢世帯が多く店舗が少ないもの）が属す地区を特定

ニュータウンポイントと行政区ポリゴンの共通部分を取れば，地区（region）情報が付与できます．

ニュータウンポイントを地区毎に分ける

変数regionの値毎に，別々のポイントデータを作成します．

地区毎に経路探索

作成したポイントデータ毎に，経路探索を行います．

探索された経路を可視化

経路を可視化します．

【補足1】テキストデータを空間データに変換

今回用いる位置情報に関するデータは，全てGeoJSON形式になっています．一方，位置情報データの中には，テキスト形式の座標値（例えば緯度経度）のみが与えられたものも多いです．上で読み込んだ店舗情報テーブルshop_tabがそれに該当します．変数lonは経度，latは緯度です．

head(shop_tab)
##      osm_id           name        shop      lon      lat
## 1 312503423 マリノアシティ        mall 130.3226 33.59501
## 2 473575921       ローソン convenience 130.4237 33.58691
## 3 477366418   ミニストップ convenience 130.3446 33.65198
## 4 485707416       BOX TOWN         yes 130.4168 33.61967
## 5 485707433 ゆめタウン博多 supermarket 130.4120 33.61187
## 6 495981699         サニー supermarket 130.4264 33.58202

座標値から空間データを作成するには，sfのst_as_sf関数を用います．引数coordsにはX,Y軸の順番で座標値変数を，crsには座標系を指定します．座標系に関する情報は，テキストデータ単体を見ても分からないので，座標がどの座標系の下で付与されたものかを，予めデータの仕様書等から確認しておく必要があります．今回の座標値は，WGS84の下で付与されたいますので，EPSGコードを4326と指定します．

shop_tab_p <- shop_tab %>%
  sf::st_as_sf(coords=c("lon","lat"),crs=sf::st_crs(4326))
head(shop_tab_p)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 130.3226 ymin: 33.58202 xmax: 130.4264 ymax: 33.65198
## Geodetic CRS:  WGS 84
##      osm_id           name        shop                  geometry
## 1 312503423 マリノアシティ        mall POINT (130.3226 33.59501)
## 2 473575921       ローソン convenience POINT (130.4237 33.58691)
## 3 477366418   ミニストップ convenience POINT (130.3446 33.65198)
## 4 485707416       BOX TOWN         yes POINT (130.4168 33.61967)
## 5 485707433 ゆめタウン博多 supermarket  POINT (130.412 33.61187)
## 6 495981699         サニー supermarket POINT (130.4264 33.58202)

こちらを更に，投影座標系（WGS84/UTM54N）に変換した場合，店舗ポイントshopとgeometry属性が一致することを確認できます．WGS84/UTM54NのEPSGコードは，32654です．

shop_tab_p_xy <- shop_tab_p %>%
  sf::st_transform(crs=sf::st_crs(32654))
head(shop_tab_p_xy)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -492920.8 ymin: 3766348 xmax: -483369.5 ymax: 3774949
## Projected CRS: WGS 84 / UTM zone 54N
##      osm_id           name        shop                  geometry
## 1 312503423 マリノアシティ        mall POINT (-492920.8 3768805)
## 2 473575921       ローソン convenience POINT (-483569.4 3766921)
## 3 477366418   ミニストップ convenience POINT (-490196.6 3774949)
## 4 485707416       BOX TOWN         yes POINT (-483832.9 3770643)
## 5 485707433 ゆめタウン博多 supermarket POINT (-484375.6 3769819)
## 6 495981699         サニー supermarket POINT (-483369.5 3766348)
head(shop)
## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -492920.8 ymin: 3766348 xmax: -483369.5 ymax: 3774949
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 3
##   osm_id    shop                   geometry
##   <chr>     <chr>               <POINT [m]>
## 1 312503423 mall        (-492920.8 3768805)
## 2 473575921 convenience (-483569.4 3766921)
## 3 477366418 convenience (-490196.6 3774949)
## 4 485707416 yes         (-483832.9 3770643)
## 5 485707433 supermarket (-484375.6 3769819)
## 6 495981699 supermarket (-483369.5 3766348)

【補足2】最寄り地物情報の付与

本編では，各メッシュから店舗や鉄道駅へのアクセス性を，メッシュ重心から1km圏内にある店舗・鉄道駅の数として定義しました．一方で，最寄り店舗や鉄道駅への距離をアクセス性の指標として用いることもできると思います．ここでは例として，各メッシュ重心から最寄りの鉄道駅と，そこまでの距離を計算する方法を示します．

メッシュ重心・鉄道駅の座標値の取得

直線距離の計算を行うためには，メッシュ重心・鉄道駅ポイントのgeometry（ポイントの座標）を，四則演算などが可能な形式の変数として格納し直す作業が必要です．変数geometryからXY座標を取り出す際には，st_coordinates関数を用います．まず，上で作成したメッシュ重心mesh_census_centを対象に実行します．

mesh_coords <- sf::st_coordinates(x=mesh_census_cent)
head(mesh_coords)
##              X       Y
## [1,] -490456.7 3753364
## [2,] -492066.2 3754941
## [3,] -491529.7 3754415
## [4,] -490945.1 3754355
## [5,] -493723.4 3756053
## [6,] -492602.5 3755467

続いて，取り出したXY座標を変数に持つメッシュ重心を新たに作成します．

mesh_census_cent_coords <- mesh_census_cent %>%
  dplyr::mutate(X_mesh=mesh_coords[,"X"],
                Y_mesh=mesh_coords[,"Y"])

ここまでの処理は，mutate関数の中でst_coordinates関数を用いることで，一括で行うことが可能です．鉄道駅ポイントについて実践します．

sta_coords <- sta %>%
  dplyr::mutate(X_sta=sf::st_coordinates(.)[,"X"],
                Y_sta=sf::st_coordinates(.)[,"Y"])

最寄り鉄道駅の探索

最寄り地物の探索とその属性の結合は，sfのst_join関数を用いて行えます．その際，引数joinにはst_nearest_featureを指定することに注意してください．

mesh_census_cent_near_sta <- mesh_census_cent_coords %>%
  sf::st_join(y=sta_coords,join=st_nearest_feature)
head(mesh_census_cent_near_sta)
## Simple feature collection with 6 features and 11 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -493723.4 ymin: 3753364 xmax: -490456.7 ymax: 3756053
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 12
##   KEY_CODE             geometry    hh eld_hh_rate   X_mesh   Y_mesh   SID
##   <chr>             <POINT [m]> <dbl>       <dbl>    <dbl>    <dbl> <int>
## 1 503012591 (-490456.7 3753364)     5       0.4   -490457. 3753364.  8960
## 2 503012674 (-492066.2 3754941)    19       0.368 -492066. 3754941.  9078
## 3 503012681 (-491529.7 3754415)    51       0.412 -491530. 3754415.  9078
## 4 503012682 (-490945.1 3754355)     5       0.4   -490945. 3754355.  9078
## 5 503012763 (-493723.4 3756053)    11       0.636 -493723. 3756053.  9078
## 6 503012771 (-492602.5 3755467)     5       0.4   -492603. 3755467.  9078
## # ℹ 5 more variables: lname <chr>, company <chr>, sname <chr>, X_sta <dbl>,
## #   Y_sta <dbl>

各メッシュ重心には最寄り駅の座標値が空間結合された状態です．座標値は投影座標系の下で付与されているので，それらをもとに最寄り駅までの直線距離が簡単に計算できます．

mesh_census_cent_near_sta_dist <- mesh_census_cent_near_sta %>%
  dplyr::mutate(d_sta=sqrt((X_sta-X_mesh)^2+(Y_sta-Y_mesh)^2))
head(mesh_census_cent_near_sta_dist)
## Simple feature collection with 6 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -493723.4 ymin: 3753364 xmax: -490456.7 ymax: 3756053
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 13
##   KEY_CODE             geometry    hh eld_hh_rate   X_mesh   Y_mesh   SID
##   <chr>             <POINT [m]> <dbl>       <dbl>    <dbl>    <dbl> <int>
## 1 503012591 (-490456.7 3753364)     5       0.4   -490457. 3753364.  8960
## 2 503012674 (-492066.2 3754941)    19       0.368 -492066. 3754941.  9078
## 3 503012681 (-491529.7 3754415)    51       0.412 -491530. 3754415.  9078
## 4 503012682 (-490945.1 3754355)     5       0.4   -490945. 3754355.  9078
## 5 503012763 (-493723.4 3756053)    11       0.636 -493723. 3756053.  9078
## 6 503012771 (-492602.5 3755467)     5       0.4   -492603. 3755467.  9078
## # ℹ 6 more variables: lname <chr>, company <chr>, sname <chr>, X_sta <dbl>,
## #   Y_sta <dbl>, d_sta <dbl>

最寄り鉄道駅距離の可視化

計算された距離変数d_staをメッシュに結合し，地図上に可視化します．そのためにまず，メッシュ重心ポイントのgeometryを落とした上で，メッシュコードと距離変数だけを残します．

mesh_census_cent_near_sta_dist_df <- mesh_census_cent_near_sta_dist %>%
  sf::st_drop_geometry() %>%
  dplyr::select(KEY_CODE,d_sta)

これを，変数KEY_CODEをキーにメッシュポリゴンに結合します．

mesh_census_dshop_ref_t_acc_dist <- mesh_census_dshop_ref_t_acc %>%
  dplyr::left_join(y=mesh_census_cent_near_sta_dist_df,by="KEY_CODE")

最後に，距離変数d_staに基づいてメッシュを色分けした可視化を行います．

mapview::mapview(x=mesh_census_dshop_ref_t_acc_dist,
                 zcol="d_sta",
                 col.regions=rev(RColorBrewer::brewer.pal(n=8,name="Spectral")),
                 color="transparent")

事例2：紛争地域の経済活動

今回の無茶振り

放送局に勤める私に，上司からこのようなメールが届きました．

昨今ウクライナ情勢が深刻さを増す中で，紛争が及ぼす経済活動への影響に関する特番を企画しました．
具体的には，この10年間において最も激しい紛争であった，シリア紛争を事例に，こうした影響を報じたいと考えています． 激戦地で撮影された写真を用いることで，経済活動への打撃を示すことを当初は考えていましたが，より客観的な統計情報を用いて影響を視覚化する必要が更にあると思います．

　そもそも，サブナショナル（1国内における地域）のレベルで，信頼性を担保しながら経済アウトカムの変遷を辿ることは，先進国を除けば平時でも難しい場合が多いです．紛争期においては，国や国際機関による統計調査の実施は尚更困難になります．
　こうした問題を乗り越えるために近年用いられているのが，夜間光データ（nighttime light data）です．夜間光は人工衛星から撮影されるため，世界の任意の場所について，同じ基準の下で経済活動を計測することができるデータです．ウクライナ紛争に際しても，NASAはTracking Night Lights in UkraineというWEBページを通じ，紛争勃発前後でのウクライナ国内における夜間光の変遷を公開しています．
　シリア紛争は，シーア派を主とするアサド政権政府軍とスンニ派の反政府軍による紛争です．イスラム国の介入等，様々な要因によって，2012年以降紛争は泥沼化し，現在も紛争は継続中です．今回は，シリア紛争勃発前（2009年）と勃発後（2018年）の2時点間で，シリア国内の夜間光量を比較し，紛争がもたらす経済活動への影響を把握します．
　夜間光にはDMSPとVIIRSの2種類があります．DMSPは1992〜2013年，VIIRSは2012年以降にデータが作成されていますが，これらデータはそれぞれ分解能等様々な側面で観測条件が異なるため，単純な前後比較が困難です．そこで今回は，以下の論文によって作成された，データ間での調整を施した夜間光データを用います．

具体的には，VIIRSのデータから，DMSPのようなデータをシミュレーションで生成したもので，例えば，2014年のロシアによるクリミア併合の経済的影響を検証した以下の論文でも用いられています．

データの読み込み

上司からのリクエストに応えるためには，シリア紛争勃発前（2009年）と勃発後（2018年）のシリア国内の夜間光データが必要そうです．また，夜間光を地域単位で集計するための，行政区域ポリゴンも必要そうです．それぞれに対応するデータとして，以下に挙げたものを用意しました．

ラスタデータ

今回用意したラスタデータは全てGeoTIFF形式です．GeoTIFFの読み込みにはterraのrast関数を用います．今回は，筆者が用意したDropboxリンクから直接ダウンロードします．

#夜間光（2009年）
ntl_09_sy <- terra::rast(x="https://www.dropbox.com/s/oekneyex34209sd/ntl_09_sy.tif?dl=1")
#夜間光（2018年）
ntl_18_sy <- terra::rast(x="https://www.dropbox.com/s/2f6td9je5daw6qz/ntl_18_sy.tif?dl=1")

読み込まれた2009年データを表示します．重要な情報は以下の通りです．

  • extent：データの空間範囲．
  • coord. ref.：データの座標系．
  • name：セルの値が入った変数の名前．
ntl_09_sy
## class       : SpatRaster 
## dimensions  : 612, 876, 1  (nrow, ncol, nlyr)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 35, 42.3, 32.3, 37.4  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : ntl_09_sy.tif?dl=1 
## varname     : ntl_09_sy 
## name        : Harmonized_DN_NTL_2009_calDMSP 
## min value   :                              0 
## max value   :                             59

セルの値は，terraのvalues関数を用いると直接表示できます．今回はセル数が膨大で，分析上直接セルの値を編集しないため，実行例は示しません．

ベクタデータ

行政区域のベクタデータを読み込みます．

#行政区域
ab_sy <- sf::read_sf(dsn="https://www.dropbox.com/s/lauux65p232v4s9/ab_sy.geojson?dl=1")

ラスタデータを可視化

読み込んだラスタデータをとりあえず可視化したい時には，terraのplot関数を用います．

#夜間光量データ（シリア・2009年）
terra::plot(ntl_09_sy)

ポリゴンによる夜間光ラスタの集計

2009・2018年のシリアの夜間光量データを，行政区域単位で合計します．各ポリゴン内に含まれる，ラスタのセルの値を集計する際には，terraのextract関数を用います．引数xには集計したいラスタを，yには集計単位となるポリゴンを（sfオブジェクトの場合，vect関数による変換が必要），funは集計方法（他にはmean, min, maxなどが利用可能），na.rmはNAのセルを除くかを指定します．

#2009年
ab_sy_ntl_09 <- terra::extract(x=ntl_09_sy,
                               y=terra::vect(ab_sy),
                               fun="sum",
                               na.rm=T) %>%
  #扱いやすい列名に変更
  dplyr::rename(ntl_09=Harmonized_DN_NTL_2009_calDMSP)
#2018年
ab_sy_ntl_18 <- terra::extract(x=ntl_18_sy,
                               y=terra::vect(ab_sy),
                               fun="sum",
                               na.rm=T) %>%
  #扱いやすい列名に変更
  dplyr::rename(ntl_18=Harmonized_DN_NTL_2018_simVIIRS)

集計された夜間光量の合計値を，シリアの行政区域ポリゴンに結合します．

ab_sy_ntl_09_18 <- ab_sy %>%
  #2009年の合計値を結合
  dplyr::left_join(y=ab_sy_ntl_09,by="ID") %>%
  #2018年の合計値を結合
  dplyr::left_join(y=ab_sy_ntl_18,by="ID")

2009・2018年の各年について，集計された夜間光量をプロットします．ポリゴンデータの可視化は，mapview関数で行えます．

#2009年
mapview::mapview(x=ab_sy_ntl_09_18,
                 zcol="ntl_09",
                 color="transparent")
#2018年
mapview::mapview(x=ab_sy_ntl_09_18,
                 zcol="ntl_18",
                 color="transparent")

2009-2018年の間の夜間光の変化率を可視化

2009・2018年の各年の夜間光量を結合した行政区域ポリゴンにおいて，両年間の夜間光量の変化率を計算します．

ab_sy_ntl_09_18_rc <- ab_sy_ntl_09_18 %>%
  dplyr::mutate(ntl_rc=(ntl_18-ntl_09)/ntl_09)

計算された変化率を，行政区域ポリゴン上に可視化します．シリア紛争の主な激戦地は，アレッポ（Aleppo）やラッカ（Raqqa），デリゾール（Deir ez-Zor）ですが，それら地域では軒並み夜間光量が減少していることがわかります．

mapview::mapview(x=ab_sy_ntl_09_18_rc,
                 zcol="ntl_rc",
                 color="transparent")

ここまでで，ようやく上司のリクエストに答えられそうな最低限の資料が作れた感じがします．

事例3：道路の建設

今回の無茶振り

建設コンサルティング会社に勤める私に，上司からこのようなメールが届きました．

昨今の中国とベトナムとの交易量増加に伴い，ベトナム北部地域と首都ハノイを結ぶ道路建設の話が持ち上がっています．
こうした道路をなるべく低い費用で建設できるルートを計画できないか，ととある機関から相談されたのですが，最大のネックは北部の険しい地形です．
なるべく勾配が少ない場所を通るようなルートで建設しないと，費用が指数関数的に増えてしまいます．
そこで，ザックリで構わないのですが，険しい地形をなるべく避けて建設できるようなルートの候補を出してもらえないでしょうか．

　例えば山岳地帯で道路建設を行う際，直線距離に近いルートを計画すると，トンネルを何個も設ける必要が出てくる場合があり，建設費用が嵩んでしまいます．故に，距離が伸びたとしても，起伏がなだらかな場所をできるだけ通すようなルートを計画した方が，時には最適です．
　こうしたルートをGISを用いて求める方法は，コストパス分析と呼ばれます．コストパス分析の基本的な考え方を示したのがFigure 2です．色が赤に近いほど傾斜が急な場所を示しますが，そうした場所をなるべく避けるような2地点間のルートを探索します．

Figure 2
Figure 2：コストパス分析の考え方（出典：Wiki.GIS.com）

データの読み込み

上司からのリクエストに応えるためには，ベトナム北部地域の地形が分かるデータが必要そうです．また，北部地域の都市の場所がわかるデータも必要そうです．それぞれに対応するデータとして，以下に挙げたものを用意しました．

ラスタデータ

標高のラスタデータを読み込みます．

el_vn <- terra::rast(x="https://www.dropbox.com/s/znei298456indc3/el_vn.tif?dl=1")

ベクタデータ

行政区域のベクタデータを読み込みます．

ab_vn_n <- sf::read_sf(dsn="https://www.dropbox.com/s/1hbvuy3prechik1/ab_vn_n.geojson?dl=1")

標高データをベトナム北部のみの領域に絞る

ラスタデータを特定の範囲だけ切り抜きたい時には，terraのcrop関数を用います．切り抜き範囲の指定は引数yで行いますが，今回は，ベトナム北部の行政区域のポリゴンをすっぽり覆う四角形（bounding boxと言います）で指定します．

el_vn_n <- terra::crop(x=el_vn,y=ext(sf::st_bbox(ab_vn_n)))
terra::plot(el_vn_n)

標高データを用いれば，傾斜度を求めることも可能です．その際には，terraのterrain関数を用います．

el_vn_n_sl <- terra::terrain(x=el_vn_n,v="slope")
terra::plot(el_vn_n_sl)

それ以外にも，地形の凹凸度（ruggedness）が同じ関数を用いて計算可能です．

el_vn_n_ru <- terra::terrain(x=el_vn_n,v="TRI")
terra::plot(el_vn_n_ru)

出発地・到着地をポイントとして与える

今回は，ベトナム最北の都市Cao Bangと，首都ハノイを結ぶコスト最小経路を探索します．故にまず，それら2つの位置をポイントデータ化する必要があります．

caoBang <- ab_vn_n %>%
  #Cao Bangのポリゴンを抽出
  dplyr::filter(NAME_2=="Cao Bang") %>%
  #ポリゴンから重心を作成
  sf::st_centroid()
haNoi <- ab_vn_n %>%
  #Ba Dinh（Ha Noiの中心地区）のポリゴンを抽出
  dplyr::filter(NAME_2=="Ba Dinh") %>%
  #ポリゴンから重心を作成
  sf::st_centroid()

movecostパッケージは，sfやterraで扱うオブジェクトには対応していません．そこで，出発地・到着地ポイントと標高データを，movecostで用いることが可能な形式に変換します．

caoBang_sp <- sf::as_Spatial(from=caoBang)
haNoi_sp <- sf::as_Spatial(from=haNoi)
el_vn_n_rl <- raster::raster(el_vn_n)

変換したファイルに基づき，movecost関数を用いたコスト最小経路探索を行います．また，探索された経路を，sfオブジェクトに変換します．

#経路探索
result_walk <- movecost::movecost(dtm=el_vn_n_rl,origin=caoBang_sp,destin=haNoi_sp,move=8)

#経路をsf形式に変換
path_walk <- sf::st_as_sf(result_walk$LCPs)

最後に，出発地・到着地とともに，探索された経路を可視化します．

mapview::mapview(x=caoBang,col.region="black") +
  mapview::mapview(x=haNoi,col.region="red") +
  mapview::mapview(x=path_walk,color="black",
                   layer.name="Route")

ここまでで，ようやく上司のリクエストに答えられそうな最低限の資料が作れた感じがします．

【補足1】カテゴリ変数のラスタデータ

ラスタデータでは，数値変数に加え，カテゴリ変数を扱うことができます．ただし，ラスタデータは文字列型の変数を持てないので，因子（factor）型でカテゴリ情報を持たせる必要があります．以下では，ベトナムの土地被覆ラスタデータ（出典：地球地図全球版）を用いて，その方法を解説します．

ラスタデータの読み込み

土地被覆のラスタデータを読み込みます．min valueとmax valueが示す通り，セルの値は1から20の間の整数値で入っています．これが土地被覆のカテゴリを表すコードです．

#土地被覆（ベトナム・2013年）
lc_vn <- terra::rast(x="https://www.dropbox.com/s/ozueb6vu9lc42h8/lc_13_vn.tif?dl=1")
lc_vn
## class       : SpatRaster 
## dimensions  : 3840, 1920, 1  (nrow, ncol, nlyr)
## resolution  : 0.004166667, 0.004166667  (x, y)
## extent      : 102, 110, 8, 24  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : lc_13_vn.tif?dl=1 
## color table : 1 
## varname     : lc_13_vn 
## name        : gm_lc_v3 
## min value   :        1 
## max value   :       20

土地被覆カテゴリの表を作成

土地被覆データの仕様ページによると，コードと土地被覆カテゴリは以下のように対応しています．

Code Class Name Code Class Name
1 Broadleaf Evergreen Forest 11 Cropland
2 Broadleaf Deciduous Forest 12 Paddy field
3 Needleleaf Evergreen Forest 13 Cropland / Other Vegetation Mosaic
4 Needleleaf Deciduous Forest 14 Mangrove
5 Mixed Forest 15 Wetland
6 Tree Open 16 Bare area,consolidated(gravel,rock)
7 Shrub 17 Bare area,unconsolidated (sand)
8 Herbaceous 18 Urban
9 Herbaceous with Sparse Tree/Shrub 19 Snow / Ice
10 Sparse vegetation 20 Water bodies

この対応表をR上で作成します．1列目にはコード，2列目には土地被覆カテゴリを変数として持ちます．

id_tab <- as.data.frame(c(1:20)) %>%
  dplyr::rename(idx=`c(1:20)`) %>%
  dplyr::mutate(nam=c("Broadleaf_Evergreen_Forest",
                      "Broadleaf_Deciduous_Forest",
                      "Needleleaf_Evergreen_Forest",
                      "Needleleaf_Deciduous_Forest",
                      "Mixed Forest",
                      "Tree_Open",
                      "Shrub",
                      "Herbaceous",
                      "Herbaceous_with_Sparse_Tree_Shrub",
                      "Sparse_vegetation",
                      "Cropland",
                      "Paddy_field",
                      "Cropland_Other_Vegetation_Mosaic",
                      "Mangrove",
                      "Wetland",
                      "Bare_area_consolidated",
                      "Bare_area_unconsolidated",
                      "Urban",
                      "Snow_Ice",
                      "Water_bodies"))
head(id_tab)
##   idx                         nam
## 1   1  Broadleaf_Evergreen_Forest
## 2   2  Broadleaf_Deciduous_Forest
## 3   3 Needleleaf_Evergreen_Forest
## 4   4 Needleleaf_Deciduous_Forest
## 5   5                Mixed Forest
## 6   6                   Tree_Open

上の対応表をラスタデータに関連づける（数値に土地被覆カテゴリをラベルとして与える）には，levels関数を用います．

levels(lc_vn) <- id_tab

この状態で，例えば行政区域毎に，extract関数を用いて土地被覆カテゴリ別のセル数をカウントすることができます．

#集計
ab_vn_n_lc <- terra::extract(x=lc_vn,
                             y=vect(ab_vn_n),
                             na.rm=T)
head(ab_vn_n_lc)
##   ID                              nam
## 1  1                      Paddy_field
## 2  1                      Paddy_field
## 3  1                      Paddy_field
## 4  1                      Paddy_field
## 5  1 Cropland_Other_Vegetation_Mosaic
## 6  1                      Paddy_field

集計結果から，行政区域毎に土地被覆シェアを計算することも可能です．コードは少々込み入っていますが，以下のような手続きを行っています．

  • 行政区域別カテゴリ別にセルの数を合計
  • 各カテゴリのセルのシェアを計算
  • 計算結果をlong型からwide型に変形
ab_vn_n_lc_sum <- ab_vn_n_lc %>%
  #行政区域（ID）と土地被覆カテゴリ変数でグループ化
  dplyr::group_by(ID,nam) %>%
  #レコード数を集計
  dplyr::summarise(count=n()) %>%
  #グループ化を解除
  dplyr::ungroup() %>%
  #行政区域（ID）でグループ化
  dplyr::group_by(ID) %>%
  #行政区域毎に，各カテゴリのセルのシェアを計算
  dplyr::mutate(prop=count/sum(count)) %>%
  #グループ化を解除
  dplyr::ungroup() %>%
  #データをlongからwideに変形
  tidyr::pivot_wider(id_cols=ID,
                     names_from=nam,
                     values_from=prop) %>%
  #NAを0に置き換え
  dplyr::mutate(across(everything(),
                       ~ifelse(!is.na(.),.,0)))
head(ab_vn_n_lc_sum)
## # A tibble: 6 × 17
##      ID Sparse_vegetation Cropland Paddy_field Cropland_Other_Vegetation…¹ Urban
##   <dbl>             <dbl>    <dbl>       <dbl>                       <dbl> <dbl>
## 1     1          0.0190    0.120        0.437                       0.0380 0.386
## 2     2          0         0.00288      0.935                       0.0615 0    
## 3     3          0.000658  0.05         0.326                       0.218  0    
## 4     4          0         0.0803       0.0918                      0.330  0    
## 5     5          0         0.0478       0.660                       0.269  0    
## 6     6          0         0.0120       0.0150                      0.119  0    
## # ℹ abbreviated name: ¹​Cropland_Other_Vegetation_Mosaic
## # ℹ 11 more variables: Tree_Open <dbl>, Broadleaf_Evergreen_Forest <dbl>,
## #   Broadleaf_Deciduous_Forest <dbl>, Needleleaf_Evergreen_Forest <dbl>,
## #   Needleleaf_Deciduous_Forest <dbl>, `Mixed Forest` <dbl>, Herbaceous <dbl>,
## #   Water_bodies <dbl>, Shrub <dbl>, Bare_area_consolidated <dbl>,
## #   Mangrove <dbl>