R で GIS データの読み込みとプロット

一昨年の目標に R で空間統計処理したいというのがありました。去年の目標は今年こそ。です。

必要なライブラリの読み込み

近年の実状を踏まえ、データ処理やプロットは tidyverse (ggplot2) を前提にし、地理空間情報の管理は sf パッケージを使うことにします。

library(tidyverse)
library(sf)

データの読み込み

UTF-8 な Shapefile

サンプルデータは Natural Earth の 1:50m Cultural Vectors countries を使用しました。

world <- st_read('./ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp')

# データの確認
head(world$NAME_JA)

[1] ジンバブエ ザンビア   イエメン   ベトナム   ベネズエラ バチカン  
241 Levels: アイスランド アイルランド ... 北マリアナ諸島

問題なく読み込まれています。

Shift_JIS な Shapefile

農林水産省が提供する2015年農林業センサスにおける農業集落境界ページで配布している福井県の市区町村境界データを使用します。

fukui.city <- st_read('./MA0003_2015_2015_18/city.shp')

# データの確認
head(fukui.city$CITY_NAME)

[1] 福井市 敦賀市 小浜市 大野市 勝山市 鯖江市
17 Levels: あわら市 おおい町 永平寺町 越前市 越前町 高浜町 坂井市 ... 福井市

なぜか Shift_JIS なデータでも正しく読み込まれています。

確実に読み込みたい場合は options を使用しエンコーディングを指定します。 st_read は内部 GDAL/OGR を使用してファイルの読み込みを行っているようなので、取り得るオプションは OGR のドライバ（ ESRI Shapefile ドライバなど）に準じます。

fukui.city <- st_read('./MA0003_2015_2015_18/city.shp', options='ENCODING=cp932')

UTF-8 な CSV

鯖江市オープンデータより消火栓データ（CSV）をダウンロード後 UTF-8 に変換した 200.utf8.csv を使用します。

hydrant <- st_read('./200.utf8.csv', options=c('X_POSSIBLE_NAMES=経度', 'Y_POSSIBLE_NAMES=緯度'), crs=st_crs(4612))

# データの確認
head(hydrant$住所)

[1] 小黒町1丁目9-27    小黒町1丁目14-14-1 小黒町1丁目14-8    小黒町3丁目12-11  
[5] 小黒町1丁目8-17    小黒町1丁目10     
2395 Levels: つつじヶ丘町1-20 つつじケ丘町11-17 ... 莇生田町34-2

CSV データを地理空間情報データとして扱うため、 options で座標として認識させる列名を指定しています。 WKT 形式の値でも可能です。詳しくは OGR の CSV ドライバ。

Shift_JIS な CSV

GDAL/OGR 自体が CSV の読み込みは UTF-8 を前提としており、エンコーディングオプションなども用意されていないため、 st_read を使って直接読み込むことはできない（？）ようです。

そこで、一旦別の関数で CSV として読み込んで、それを sf データにします。

hydrant <- read_csv('200.csv', locale=locale(encoding='cp932')) %>%
    st_as_sf(coords=c('経度', '緯度'), crs=st_crs(4612))

# データの確認
head(hydrant$住所)

[1] "小黒町1丁目9-27"    "小黒町1丁目14-14-1" "小黒町1丁目14-8"   
[4] "小黒町3丁目12-11"   "小黒町1丁目8-17"    "小黒町1丁目10"

地図のプロット

基本形

ggplot() +
    geom_sf(data=world)

プロットする地図の投影法

ggplot() +
    geom_sf(data=world) +
    coord_sf(crs=st_crs(3857))

プロットする地図の範囲

メルカトル図法は極圏を表示しようとすると縦長になってしまいますので、 xlim, ylim で表示範囲を指定します。ただし、投影後の座標で指定する必要があります。

ここでは経緯度から投影後の座標を計算することにします。

xlim <- (st_point(c(180,0)) %>%
    st_sfc(crs=4326) %>%
    st_transform(st_crs(3857)) %>%
    as_vector())[1] %>%
    abs()

EPSG:4326 で (180, 0) の座標を EPSG:3857 に投影変換し、 X 座標を取得しています。なぜか負数になるのですが、理由がいまいち不明なので abs() を使用して確実に正数にしています。

ylim についても計算してもよいですが、ここでは正方形にすることとします。

ylim <- xlim

ggplot() +
    geom_sf(data=world) +
    coord_sf(crs=st_crs(3857), xlim=c(-xlim,xlim), ylim=c(-ylim,ylim))

複数データの読み込み

ggplot() +
    geom_sf(data=fukui.city %>% filter(`CITY_NAME`=='鯖江市')) +
    geom_sf(data=hydrant)

プロットする地図の投影法を変え、目盛の単位も変える

ggplot() +
    geom_sf(data=fukui.city %>% filter(`CITY_NAME`=='鯖江市')) +
    geom_sf(data=hydrant) +
    coord_sf(crs=st_crs(6674), datum=st_crs(6674))

目盛軸の間隔を変える

ggplot() +
    geom_sf(data=fukui.city %>% filter(`CITY_NAME`=='鯖江市')) +
    geom_sf(data=hydrant) +
    coord_sf(crs=st_crs(6674), datum=st_crs(6674)) +
    scale_x_continuous(breaks=seq(10000, 30000, 2000)) +
    scale_y_continuous(breaks=seq(-8000,0,2000))