library(dplyr) # データ操作
library(sf) # 地理空間データの操作
library(mapview) # 地図のインタラクティブな操作
library(ggplot2) # データ可視化
library(rnaturalearth) # パブリックドメインで利用可能な県単位のポリゴンデータ
library(rnaturalearthhires)地図へのマッピング
多くの統計データには、都道府県名やメッシュ、緯度経度の座標といった地理的な空間情報を含みます。 このような地理空間情報をもつデータ（地理空間データ）は地図上に投影することで内容の理解を深めることが期待できます。 また、地理空間データの操作により分析の幅を広げることができます。 Rの地理空間データを使って、これらの処理を行いましょう。
1. パッケージの読み込み
ここでもいくつかのパッケージを使います。 関数を利用するパッケージを次のコマンドで読み込みます。
2. 地理空間データの用意
オープンデータの活用の中で取り上げた「緊急避難場所 (徳島県)」を再び利用します。 このデータは県内の市町村が指定する洪水災害発生時に利用可能な施設に関しての位置情報を記録したものでした。 このデータに含まれる位置情報（緯度、経度）をもとに、Rで地理空間データを作成してみます。
sfパッケージのst_as_sf()関数を使って、データフレーム中の緯度 latitudeと経度 longitudeの列の値をもとに地理空間情報データ（ポイント）に変換します。
# 緯度経度の座標からポイントデータを生成
sf_hinanjyo <-
st_as_sf(df_hinanjyo,
coords = c("longitude", "latitude"),
crs = 4326)
glimpse(sf_hinanjyo)Rows: 701
Columns: 5
$ city <chr> "板野郡板野町", "阿南市", "阿南市", "阿南市", "阿南市", "阿南…
$ title <chr> "板野町防災ステーション", "新野小学校", "新野東小学校", "阿南…
$ address <chr> "徳島県板野郡板野町川端字新手崎18-1", "阿南市新野町南宮ノ久保…
$ type <chr> "指定緊急避難場所（法指定）", "指定緊急避難場所（法指定）", "…
$ geometry <POINT [°]> POINT (134.4729 34.13761), POINT (134.5811 33.84699), P…geometry列に各施設の位置がポイントデータとして記録されています。
続いて、徳島県の形状を示すためのポリゴンデータを準備します。 ここではパブリックドメインで使用できるNatural EarthのデータをRから取得するrnaturalearthパッケージを用います。 次のコマンドで四国４県のポリゴンデータを用意します。
# rnaturalearthから四国のポリゴンを準備する
ne_jpn_shikoku <-
rnaturalearth::ne_states(country = "Japan", returnclass = "sf") |>
filter(region == "Shikoku") |>
select(iso_3166_2, name)
ne_jpn_shikoku# A tibble: 4 × 3
iso_3166_2 name geometry
<chr> <chr> <MULTIPOLYGON [°]>
1 JP-36 Tokushima (((134.4424 34.20827, 134.4686 34.20746, 134.4952 34.224…
2 JP-37 Kagawa (((133.5919 34.02381, 133.6209 34.04556, 133.634 34.0607…
3 JP-38 Ehime (((132.6399 32.90892, 132.6272 32.90632, 132.6184 32.906…
4 JP-39 Kochi (((134.2957 33.53032, 134.2732 33.508, 134.2446 33.46075…このほか、行政区域のデータには国土数値情報ダウンロードサービスの行政区域データなどが利用できます。国土数値情報の行政区域データは市区町村別にポリゴンデータがあたえられているため、細かな地域の可視化や分析に役立ちます。
3. 地図の作成
Rでの地理空間データの可視化表現として、インタラクティブに操作できる地図、グラフの作成の中で扱ったggplot2による静的な図の作成ができます。
インタラクティブな地図の操作
避難場所の位置を地図上で確認します。
mapviewパッケージはインタラクティブに操作できる地図機能を提供します。 mapview()関数に対象の地理空間データを与えて実行すると、自在に動かせる地図の画面が表示されます。 以下に出力される地図を操作（移動、拡大・縮小、レイヤの変更、アイコンのクリック）をしてみましょう。
mapview(sf_hinanjyo)# 市町村別の塗り分け
# mapview(sf_hinanjyo, zcol = "city")同様の出力として、避難所等の場所を地図上にマッピングするアプリケーションを徳島県が公開しています。関心のある方はRでの出力と県のページを比較してみてください。
静的な地図の描画とレイヤの重ね合わせ
静的な地図はggplot2パッケージのgeom_sf()関数により生成可能です。 これにより、徳島県の形状を表すポリゴンデータ、避難場所の位置を示すポイントデータの2つを重ね合わせた地図を作成できます。
ggplot() +
geom_sf(data = filter(ne_jpn_shikoku, name == "Tokushima")) +
geom_sf(data = sf_hinanjyo,
aes(color = city),
show.legend = FALSE) +
coord_sf() +
labs(title = "徳島県緊急避難場所の位置")
4. 地理空間データの処理
地理空間データの処理はsfパッケージを介して行います。 例えば、ポリゴンデータから面積を計算する、2地点（ポイント）間の距離を求める、などが可能です。
# 面積の算出
st_area(ne_jpn_shikoku)Units: [m^2]
[1] 4249555713 1853829862 5633714966 7091706264# 単位の変換
units::set_units(st_area(ne_jpn_shikoku), km^2)Units: [km^2]
[1] 4249.556 1853.830 5633.715 7091.706# 距離の計算
st_distance(
sf_hinanjyo[1, ], # 板野町防災ステーション (徳島県板野郡板野町川端字新手崎18-1)
sf_hinanjyo[10, ] # 阿南市クリーンピュア (阿南市熊谷町定方44)
)Units: [m]
[,1]
[1,] 26405.25# ポリゴンの重心点
st_centroid(ne_jpn_shikoku)Warning in st_centroid.sf(ne_jpn_shikoku): st_centroid assumes attributes are
constant over geometries of x# A tibble: 4 × 3
iso_3166_2 name geometry
<chr> <chr> <POINT [°]>
1 JP-36 Tokushima (134.2332 33.91233)
2 JP-37 Kagawa (134.0072 34.23647)
3 JP-38 Ehime (132.8498 33.61348)
4 JP-39 Kochi (133.3479 33.40644)sfパッケージで求めた値は、dplyrパッケージのmutate()関数と組み合わせることでデータフレームの列として格納できます。
mutate(ne_jpn_shikoku,
area = units::set_units(st_area(ne_jpn_shikoku), km^2))# A tibble: 4 × 4
iso_3166_2 name geometry area
<chr> <chr> <MULTIPOLYGON [°]> [km^…
1 JP-36 Tokushima (((134.4424 34.20827, 134.4686 34.20746, 134.4952 … 4250.
2 JP-37 Kagawa (((133.5919 34.02381, 133.6209 34.04556, 133.634 3… 1854.
3 JP-38 Ehime (((132.6399 32.90892, 132.6272 32.90632, 132.6184 … 5634.
4 JP-39 Kochi (((134.2957 33.53032, 134.2732 33.508, 134.2446 33… 7092.応用例として、徳島県立総合教育センター(板野郡板野町犬伏字東谷1-7)から半径1km圏内にある避難場所を特定します。 まずは徳島県立総合教育センターの座標をst_point()関数内で定義し、ポイントデータを作成します。
x <-
# 徳島県立総合教育センターの位置
st_sfc(st_point(c(134.452485, 34.150387)),
# 座標参照系の指定
crs = 4326)
# 位置を確認
mapview(x, map.types = "OpenStreetMap")次にst_buffer()関数でポイントからバッファを生成します。 ここで引数distにバッファの大きさを与えて実行します。 今回は半径1kmのバッファとしたいので、units::set_units(1, km)を与えます。
x_buffer1km <-
st_buffer(x,
dist = units::set_units(1, km))
mapview(x_buffer1km, map.types = "OpenStreetMap")sf_hinanjyo_neighborhood <-
st_join(
sf_hinanjyo,
st_sf(x_buffer1km),
left = FALSE)
sf_hinanjyo_neighborhood# A tibble: 5 × 5
city title address type geometry
<chr> <chr> <chr> <chr> <POINT [°]>
1 板野町 穂波園指定通所介護事業所 板野町吹田字… 指定… (134.4555 34.1458)
2 板野町 文化の館 板野町犬伏字… 指定… (134.4534 34.14513)
3 板野町 板野町体育センター 板野町吹田字… 指定… (134.4555 34.14728)
4 板野町 犬伏老人憩の家 板野町犬伏字… 指定… (134.4506 34.14246)
5 板野町 吹田老人憩の家 板野町吹田字… 指定… (134.4608 34.14848)5地点の避難場所が示されました。 これらの避難場所と徳島県立総合教育センターの距離をst_dist()関数で求めます。
sf_hinanjyo_neighborhood <-
mutate(sf_hinanjyo_neighborhood,
distance = st_distance(geometry, x, by_element = FALSE))最後に、避難場所の位置関係も確認しておきましょう。 mapviewパッケージではggplot2パッケージのように描画対象のオブジェクトを+演算子を使って レイヤとして重ねられます。
mapview(x_buffer1km) +
mapview(sf_hinanjyo_neighborhood,
zcol = "distance")大きなバッファの内側に対象の5地点の避難場所が含まれることが確認できました。
5. まとめ
演習
任意の地点の座標（緯度経度）を調べ、その地点から半径1km圏内にある避難場所を調べてみよう。