シェープファイルを読む

[2020-09-17] Rが4.xになったことと，国土数値情報の令和2年（2020年）版が出たことと，WindowsでRをお使いの方からフィードバックをいただけたことにより，このページと集約geojsonデータを改訂しました。

はじめに

シェープファイル（Shapefile）は，地図を表すベクトルデータ形式である。Esri が提唱した。拡張子 shp，shx，dbf の三つのファイルから成る。ほかにも投影法についての情報を持つ prj や，メタデータの xml などを含むことがある。

日本の行政区域を表すシェープファイルは，国土数値情報サイトの「行政区域」や，e-Stat の「地図で見る」→「境界データダウンロード」→「小地域」→「国勢調査」→「小地域」→「世界測地系緯度経度・Shape形式」などで入手できる。

国土数値のファイルをダウンロードする際に kokudosuuchi パッケージが便利である。参考：Rから国土数値情報ダウンロードサービス Web APIを使うパッケージkokudosuuchiをCRANで公開しました

シェープファイルを読むRのパッケージはいろいろあるが，ここでは sf（simple features）パッケージを用いる。

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

シェープファイルを読む関数はいくつかあるが read_sf() が良さそうである。

三重県のシェープファイル

ここでは国土数値情報の「行政区域」の三重県を使ってみる。ファイル名は N03-200101_24_GML.zip である。これをカレントディレクトリの N03-200101_24_GML に展開し，Rでシェープファイルを読み込む。中身がシフトJIS（CP932）のため，エンコーディングの指定が必要である。ファイル名を指定する場合は shp のファイル名だけでよいし，shp，shx，dbf の入っているディレクトリを指定するだけでもよい（カレントディレクトリなら "."）。ただし，WindowsのRでファイル名を指定しないと動かなかったというご報告をいただいたので，以下ではshpファイル名を指定するようにした。

mie = read_sf("N03-200101_24_GML/N03-20_24_200101.shp", options="ENCODING=CP932")

中を見てみよう：

mie
Simple feature collection with 7480 features and 5 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 135.8532 ymin: 33.72286 xmax: 136.9901 ymax: 35.25765
geographic CRS: JGD2011
# A tibble: 7,480 x 6
   N03_001 N03_002 N03_003 N03_004  N03_007                             geometry
   <chr>   <chr>   <chr>   <chr>    <chr>                          <POLYGON [°]>
 1 三重県  NA      NA      津市     24201   ((136.2062 34.44868, 136.2065 34.44…
 2 三重県  NA      NA      津市     24201   ((136.5276 34.7106, 136.5277 34.710…
 3 三重県  NA      NA      四日市市 24202   ((136.625 34.90368, 136.6249 34.903…
 4 三重県  NA      NA      四日市市 24202   ((136.6458 34.93301, 136.6458 34.93…
 5 三重県  NA      NA      四日市市 24202   ((136.6466 34.94489, 136.6466 34.94…
 6 三重県  NA      NA      四日市市 24202   ((136.6583 34.95221, 136.6573 34.95…
 7 三重県  NA      NA      四日市市 24202   ((136.6522 34.98956, 136.6522 34.98…
 8 三重県  NA      NA      四日市市 24202   ((136.6527 34.994, 136.6528 34.9940…
 9 三重県  NA      NA      四日市市 24202   ((136.6833 34.98159, 136.6831 34.98…
10 三重県  NA      NA      四日市市 24202   ((136.6742 35, 136.673 34.99794, 13…
# … with 7,470 more rows
table(mie$N03_003)

  員弁郡   桑名郡   三重郡   多気郡   度会郡 南牟婁郡 北牟婁郡 
       1        1        6        7     2301       27     1362

N03_007 は全国地方公共団体コードの市区町村コードからチェックディジットを除いたものである。

もし読み込んでから文字化けがわかった場合，次のようにしてUTF-8に変換できる：

mie$N03_001 = iconv(mie$N03_001, from="CP932", to="UTF-8")
mie$N03_003 = iconv(mie$N03_003, from="CP932", to="UTF-8")
mie$N03_004 = iconv(mie$N03_004, from="CP932", to="UTF-8")

これをプロットするために

plot(mie)

と打ち込むと，N03_001，N03_002，N03_003，N03_004，N03_007 のそれぞれについて三重県の地図が出力される。N03_001 はすべて「三重県」であるのですべて同じ色で，N03_002 はすべて NA であるので白地図で，N03_003 は郡部だけ色分けして，N03_004 と N03_007 は行政区分ごとに色分けして出力される。

ポリゴンデータの部分は mie$geometry でもアクセスできるが，より一般的には st_geometry(mie) とする：

st_geometry(mie)
Geometry set for 7480 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 135.8532 ymin: 33.72286 xmax: 136.9901 ymax: 35.25765
geographic CRS: JGD2011
First 5 geometries:
POLYGON ((136.2062 34.44868, 136.2065 34.44904,...
POLYGON ((136.5276 34.7106, 136.5277 34.71062, ...
POLYGON ((136.625 34.90368, 136.6249 34.90368, ...
POLYGON ((136.6458 34.93301, 136.6458 34.93301,...
POLYGON ((136.6466 34.94489, 136.6466 34.94488,...

したがって，

plot(st_geometry(mie))

とだけ打ち込めば行政区分の白地図ができる。余白を狭くしたければあらかじめ par(mar=c(0,0,0,0)) と打ち込んでおけばよい。また，区分なしの全体を描きたければ，次のようにする：

plot(st_union(mie))

市町村単位の情報を描きこむには，1行が1市町村に対応するようにデータを集約するほうが便利である：

mie2 = aggregate(mie, list(mie$N03_007), unique)

こうすれば，jpndistrictによる地図で描いたような図が描きやすい：

x = c(11,3,100,6,34,25,1,4,40,91,0,100,0,84,100,100,71,58,29,18,38,100,100,100,83,100,10,0,0)
cols = colorRamp(c("#66ccff", "white", "#ff9900"))
plot(st_geometry(mie2), col=rgb(cols(x/100)/255))

各市町村の重心座標を求め，そこに市町村名を書く：

c = st_centroid(st_geometry(mie2))
text(st_coordinates(c), mie2$N03_004, cex=0.5)

逆に，最初読み込んだままの行を生かすことも可能である：

t = unique(data.frame(code=mie$N03_007, name=mie$N03_004))
t$x = x
mie3 = merge(mie, t, by.x="N03_007", by.y="code")
plot(st_geometry(mie3), col=rgb(cols(mie3$x/100)/255))

全国のシェープファイル

国土数値情報の N03-200101_GML.zip をいただいてきて展開する。こちらのほうはシフトJISではないのでエンコーディングのオプションは不要であった。

japan = read_sf("N03-200101_GML/N03-20_200101.shp")

これをコードの上2桁で集約すると，1行が1都道府県になる（とても時間がかかる）：

japan2 = aggregate(japan, list(substr(japan$N03_007, 1, 2)), head, n=1)

これをそのままプロットすると非常に重く，しかも境界が複雑で見栄えが悪い。そこで境界を単純化する。そのためには sf パッケージの st_simplify() が使えるが，より良いアルゴリズムが rmapshaper パッケージの ms_simplify() に実装されている。データ量を約 1/1000 に縮めてみよう：

install.packages("rmapshaper")
library(rmapshaper)
japan3 = ms_simplify(japan2, keep=0.001, keep_shapes=TRUE)

sf パッケージの st_simplify() でも試してみよう。こちらは緯度経度データを単純化しようとすると警告が出る。無視してもよいが，km 単位に直してから 1 km 精度で単純化してみる：

japan2s = st_transform(japan2, "+proj=utm +zone=54 +datum=WGS84 +units=km")
japan3s = st_simplify(japan2s, dTolerance=1)
japan3s = st_transform(japan3s, "+proj=longlat +ellps=GRS80")

余計な列を省き，列の名前を付け替える：

japan = japan3[ ,1:2]
names(japan) = c("code","pref","geometry")

できたデータをシェープファイル形式で保存するには

st_write(japan, "japan.shp", layer_options="ENCODING=UTF-8")

とすればよいが，ここではGeoJSON形式で保存する。ファイルが一つにまとまっておりネット越しに使うのも楽である：

st_write(japan, "japan.geojson")

できた japan.geojson は公開しておく。次のように読む：

japan = read_sf("https://okumuralab.org/~okumura/stat/data/japan.geojson")

確認には plot(st_geometry(japan)) としてみればよいが，ここでは投影法の確認も兼ねて，緯度・経度の線（graticule）や軸（axes）も描きこんでみよう：

plot(st_geometry(japan), graticule=TRUE, axes=TRUE)

座標系をUTM（Universal Transverse Mercator）に変えてみよう：

j = st_transform(japan, "+proj=utm +zone=54 +datum=WGS84 +units=km") # 52:九州 53:西日本 54:東日本
plot(st_geometry(j), graticule=TRUE, axes=TRUE)

例として日本地図: 普通教室のエアコン設置率と同じことをしてみよう：

x = read.csv("https://okumuralab.org/~okumura/stat/data/aircon_shochu.csv", fileEncoding="UTF-8")
cols = colorRamp(c("#ff2800","white","#0041ff"))
pal = rgb(cols((0:101)/101)/255)
j = st_transform(japan, "+proj=utm +zone=54 +datum=WGS84 +units=km") # 52:九州 53:西日本 54:東日本
j$aircon = x[,4]   # 範囲 [0,100]
par(mar=c(0,0,0,0))
par(mgp=c(2,0.8,0))
plot(j["aircon"], lwd=0.3, pal=pal, breaks=-1:101, nbreaks=103, key.length=0.4, key.pos=4, main="")

もうちょっと単純化して，沖縄県はインセットに入れてみる：

j = ms_simplify(j, keep=0.1, keep_shapes=TRUE)
j$geometry[[47]] = j$geometry[[47]] + c(700,1500)
plot(j["aircon"], lwd=0.3, pal=pal, breaks=-1:101, nbreaks=103, key.length=0.4, key.pos=4,
     main="", xlim=c(-600, 1050), ylim=c(3500, 5030), reset=FALSE)
lines(c(-600,-300,100,100), c(4150,4150,4450,4750))

ついでに全国を市町村レベルで集約して，約 1/100 に単純化したものを用意しておこう。重い作業なので，いったんRを再起動して最初からやり直す：

library(sf)
library(rmapshaper)
# options(stringsAsFactors=FALSE) # R 4.x ではデフォルト
japan = read_sf("N03-200101_GML/N03-20_200101.shp")
japan2 = aggregate(japan, list(japan$N03_007), head, n=1)
japan3 = ms_simplify(japan2, keep=0.01, keep_shapes=TRUE)
#
# Fatal error in , line 0
# API fatal error handler returned after process out of memory
#

今回はメモリ16GのMac miniでやったのだが，やはり R が落ちてしまう。そこで，北海道，本州，残りに分割する：

library(sf)
library(rmapshaper)
# options(stringsAsFactors=FALSE)
japan = read_sf("N03-200101_GML/N03-20_200101.shp")
japan2 = aggregate(japan, list(japan$N03_007), head, n=1)
japan2$Group.1 = NULL
japan2a = japan2[as.numeric(substr(japan2$N03_007, 1, 2)) %in% 1,]
japan3a = ms_simplify(japan2a, keep=0.01, keep_shapes=TRUE)
japan2b = japan2[as.numeric(substr(japan2$N03_007, 1, 2)) %in% 2:35,]
japan3b = ms_simplify(japan2b, keep=0.01, keep_shapes=TRUE)
japan2c = japan2[as.numeric(substr(japan2$N03_007, 1, 2)) %in% 36:47,]
japan3c = ms_simplify(japan2c, keep=0.01, keep_shapes=TRUE)
japan3 = rbind(japan3a,japan3b,japan3c)
st_write(japan3, "japan2.geojson")

できた japan2.geojson は公開しておく。次のように読む：

japan2 = read_sf("https://okumuralab.org/~okumura/stat/data/japan2.geojson", stringsAsFactors=FALSE)

これを使って三重県の白地図を描こう：

mie = subset(japan2, N03_001=="三重県")
par(mar=c(0,0,0,0))
plot(st_geometry(mie), lwd=0.3)
c = st_centroid(st_geometry(mie))
text(st_coordinates(c), mie$N03_004)

RgoogleMapsを使った放射線地図もこれを使って描くことができる：

library(sf)
library(RColorBrewer)
japan = st_read("https://okumuralab.org/~okumura/stat/data/japan2.geojson", stringsAsFactors=FALSE)
fukushima = read.csv("https://okumuralab.org/~okumura/stat/data/fukushima.csv", as.is=TRUE)
t = as.POSIXct(fukushima$datetime)
o = order(t)
cols=c("white",brewer.pal(9, "YlOrRd"),rep("black",100))
par(mar=c(0,0,0,0))
plot(st_geometry(japan2), xlim=c(139.2015,141.0421), ylim=c(36.7540,38.01304), lwd=0.3, col="#c8c8cb")
c = st_centroid(st_geometry(japan2))
points(fukushima$lon[o], fukushima$lat[o], pch=16, col=cols[floor(fukushima$radiation[o]*2)+1])
points(141.032339, 37.422778, pch="×", cex=2)
text(st_coordinates(c), japan2$N03_004, cex=0.5)

ついでに，福島第一原発と三重大学の距離（大円距離）は500kmくらいである：

f1 = st_sfc(st_point(c(141.032339,37.422778)), crs="+proj=longlat")
mieu = st_sfc(st_point(c(136.5248,34.7468)), crs="+proj=longlat")
st_distance(f1, mieu, by_element=TRUE)
502871.3 [m]