シェープファイルを入手

シェープファイルは国土数値情報　行政区域データの詳細から全国分を取得。

シェープファイルを読み込む

ダウンロードしたZIPを展開したフォルダが作業ディレクトリ内においてあって、N03-160101_GML/N03-16_160101.shpにシェープファイルがあるものと仮定する。

シェープファイルはsfパッケージを使用して読み込む。

library(sf)
# 日本地図のシェープファイルを読み込む
f <- "N03-160101_GML/N03-16_160101.shp"
japan <- st_read(f)

中身を確認するとN03_001に都道府県名が入っていることが確認できる。

> print(japan, n = 3)
Simple feature collection with 115866 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
First 3 features:
  N03_001    N03_002 N03_003 N03_004 N03_007                       geometry
1  北海道 石狩振興局  札幌市  中央区   01101 POLYGON((141.342326939 43.0...
2  北海道 石狩振興局  札幌市    北区   01102 POLYGON((141.408387224 43.1...
3  北海道 石狩振興局  札幌市    東区   01103 POLYGON((141.447070272 43.1...

都道府県ごとにポリゴンをまとめる

st_union()で一つに結合することもできるが、st_combine()で単にまとめてMULTIPOLYGONとしてしまったほうが早い。境界は残るが、実用上問題はないだろうと思う。

# 都道府県名一覧の取得
library(dplyr)
japan_p <- select(japan, N03_001)
pref <- unique(japan$N03_001)
# 都道府県ごとにポリゴンをまとめる
pref_polygons <- tapply(japan_p$geometry,
                        japan_p$N03_001, # 都道府県名でグループ化
                        st_combine)      # st_unionより早い

ポイントがどの県に含まれているのかを調べる

オブジェクトの包含関係はst_contains()で調べることができる。
st_contains(A,B)でAにBが含まれると数値1を含むリストが返ってきて、含まれないと長さ0の数値ベクトルを含むリストが帰ってくる。返り値が妙なので使い方を間違えているのかもしれない。

pref <- unique(japan$N03_001) # 都道府県名の取得
find_pref <- function(point){ # 力技で検索する
  for(p in pref){
    if(length(st_contains(pref_polygons[[p]], point)[[1]])){
      return(p)
    }
  }
  return(NA)
}

無理矢理検索しているので若干時間がかかる。

> p1 <- c(138.5, 37.1) # 新潟県
> p2 <- c(138.5, 36.8) # 長野県
> p3 <- c(138.5, 37.5) # うみのうえ
> find_pref(st_point(p1))
[1] "新潟県"
> find_pref(st_point(p2))
[1] "長野県"
> find_pref(st_point(p3))
[1] NA

ここまでやったところで

別に包含関係なんか調べなくても重心が一番近い行政区域をひけば逆ジオコーディングに近い何かになるのでは？と思った。変な形の場所があると問題が出そうだけど。

jpn_cent <- st_centroid(japan) # セントロイド計算
cent_mat <- t(matrix(unlist(jpn_cent$geometry), nrow=2)) # 地点をmatrixに

mydist <- function(p1, p2){ # 2点間のユークリッド距離
  sqrt(sum((p1 - p2)^2))
}
find_admin <- function(point){ # 最寄りの行政区域
  jpn_cent[which.min(apply(cent_mat, 1, function(x) mydist(x, point))),]
}

> find_admin(p1)
Simple feature collection with 1 feature and 5 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 138.3375 ymin: 37.12011 xmax: 138.3375 ymax: 37.12011
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=GRS80 +no_defs
      N03_001 N03_002 N03_003 N03_004 N03_007                       geometry
40776  新潟県    <NA>    <NA>  上越市   15222 POINT(138.337517708344 37.1...
> find_admin(p2)
Simple feature collection with 1 feature and 5 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 138.5106 ymin: 36.75132 xmax: 138.5106 ymax: 36.75132
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=GRS80 +no_defs
      N03_001 N03_002  N03_003  N03_004 N03_007                       geometry
47022  長野県    <NA> 下高井郡 山ノ内町   20561 POINT(138.510638855801 36.7...
> find_admin(p3)
Simple feature collection with 1 feature and 5 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 138.6165 ymin: 37.47928 xmax: 138.6165 ymax: 37.47928
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=GRS80 +no_defs
      N03_001 N03_002 N03_003 N03_004 N03_007                       geometry
39669  新潟県    <NA>    <NA>  柏崎市   15205 POINT(138.61649071978 37.47...

細かい情報が得られる上に5~10倍くらい早い。

> system.time(find_pref(st_point(p1)))
   ユーザ   システム       経過  
     1.784      0.199      1.990 
> system.time(find_pref(st_point(p2)))
   ユーザ   システム       経過  
     2.366      0.245      2.620 
> system.time(find_pref(st_point(p3)))
   ユーザ   システム       経過  
     4.437      0.466      4.921 
> system.time(find_admin(p1))
   ユーザ   システム       経過  
     0.459      0.004      0.464 
> system.time(find_admin(p2))
   ユーザ   システム       経過  
     0.437      0.002      0.441 
> system.time(find_admin(p3))
   ユーザ   システム       経過  
     0.456      0.004      0.462

追記

#mydist <- function(p1, p2){ # 2点間のユークリッド距離
#  sqrt(sum((p1 - p2)^2))
#}
find_admin <- function(point){ # 最寄りの行政区域
  #jpn_cent[which.min(apply(cent_mat, 1, function(x) mydist(x, point))),]
  jpn_cent[which.min(sqrt(rowSums(t(t(cent_mat) - point)^2))),]
}

100倍早くなった。ベクトル化大切。

> system.time(find_pref(p1))
   ユーザ   システム       経過  
     0.003      0.001      0.004