シェープファイルを読む

はじめに

シェープファイル(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)

シェープファイルを読む関数は st_read() または read_sf() である。前者はメッセージを出力し,後者は出力しない。前者はデフォルトで文字列をファクターに変換し,後者は変換しない。いずれにしても options(stringsAsFactors=FALSE) を実行しておけばファクターに変換されない。

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

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

options(stringsAsFactors=FALSE)
mie = st_read("N03-180101_24_GML/", options="ENCODING=CP932")
options:        ENCODING=CP932 
Reading layer `N03-18_24_180101' from data source `/.../.../N03-180101_24_GML' using driver `ESRI Shapefile'
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
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=GRS80 +no_defs

中を見てみよう:

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
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=GRS80 +no_defs
# 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.5282 34.71053, 136.5276 34.7…
 2 三重県  NA      NA      津市     24201   ((136.4232 34.84402, 136.4233 34.8…
 3 三重県  NA      NA      四日市市 24202   ((136.6458 34.93301, 136.6458 34.9…
 4 三重県  NA      NA      四日市市 24202   ((136.6466 34.94488, 136.6465 34.9…
 5 三重県  NA      NA      四日市市 24202   ((136.6589 34.95153, 136.6588 34.9…
 6 三重県  NA      NA      四日市市 24202   ((136.6771 34.97827, 136.6763 34.9…
 7 三重県  NA      NA      四日市市 24202   ((136.652 34.98822, 136.6519 34.98…
 8 三重県  NA      NA      四日市市 24202   ((136.6531 34.99381, 136.653 34.99…
 9 三重県  NA      NA      四日市市 24202   ((136.673 34.99794, 136.6728 34.99…
10 三重県  NA      NA      四日市市 24202   ((136.5404 35.06994, 136.5407 35.0…
# ... 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
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=GRS80 +no_defs
First 5 geometries:
POLYGON ((136.5282 34.71053, 136.5276 34.71027,...
POLYGON ((136.4232 34.84402, 136.4233 34.84401,...
POLYGON ((136.6458 34.93301, 136.6458 34.93301,...
POLYGON ((136.6466 34.94488, 136.6465 34.94479,...
POLYGON ((136.6589 34.95153, 136.6588 34.9515, ...

したがって,

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-180101_GML.zip をいただいてきて展開する。こちらのほうはシフトJISではないのでエンコーディングのオプションは不要であった。

options(stringsAsFactors=FALSE)
japan = st_read("N03-180101_GML/")
Reading layer `N03-18_180101' from data source `/.../N03-180101_GML' using driver `ESRI Shapefile'
Simple feature collection with 116929 features and 5 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 122.9337 ymin: 20.42275 xmax: 153.9868 ymax: 45.55724
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=GRS80 +no_defs

これをコードの上2桁で集約すると,1行が1都道府県になる(とても時間がかかる):

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

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

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 = st_read("https://oku.edu.mie-u.ac.jp/~okumura/stat/data/japan.geojson", stringsAsFactors=FALSE)

確認には 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)
緯度経度 UTM 54

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

x = read.csv("https://oku.edu.mie-u.ac.jp/~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 に単純化したものを用意しておこう:

library(sf)
library(rmapshaper)
options(stringsAsFactors=FALSE)
japan = st_read("N03-180101_GML/")
japan2 = aggregate(japan, list(japan$N03_007), head, n=1)
japan3 = ms_simplify(japan2, keep=0.01, keep_shapes=TRUE)
#
# Fatal error in CALL_AND_RETRY_0
# Allocation failed - process out of memory
#

データ量が多いと,このように R が落ちてしまう。そこで,北海道,本州,残りに分割する:

library(sf)
library(rmapshaper)
options(stringsAsFactors=FALSE)
japan = st_read("N03-180101_GML/")
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 = st_read("https://oku.edu.mie-u.ac.jp/~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://oku.edu.mie-u.ac.jp/~okumura/stat/data/japan2.geojson", stringsAsFactors=FALSE)
fukushima = read.csv("https://oku.edu.mie-u.ac.jp/~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)
放射線地図

参考