Rで緯度経度から都道府県・市区町村を求める

@nozma さんが書いてくれた記事の俺俺改造記事です。nozmaさんの記事では「都道府県ごとにポリゴンをまとめる」という工程を行なっていますが、それを使わない方法、日本全国のポリゴンを直接利用する方法を紹介します。

なお国土数値情報の提供する行政区域データをRから呼び出すためのjpndistrictというパッケージがありますが、これは都道府県レベルのShapefileにしか対応していないので、全国shapefileは直接ダウンロードすることになります。今回は「2014年版 行政区域データ 全国のファイル」を使います。

余談ですが、jpndistrictパッケージはCRANに登録したものの、諸事情によりアーカイブされています。この問題を解決した次のバージョンではこの記事で取り上げる都道府県・市区町村判定の関数も取り込む予定です。

library(magrittr)
library(dplyr)
library(kokudosuuchi)

## このサービスは、「国土交通省　国土数値情報（カテゴリ名）」をもとに加工者が作成
## 以下の国土数値情報ダウンロードサービスの利用約款をご確認の上ご利用ください：
## 
## http://nlftp.mlit.go.jp/ksj/other/yakkan.html

# 2014年版 行政区域データ 全国のファイルをダウンロードして解凍
dl.url <- kokudosuuchi::getKSJURL(identifier = "N03") %>%
  dplyr::filter(year == 2014, areaCode == 0) %>% use_series(zipFileUrl)

download.file(dl.url,
              destfile = "N03-150101_GML.zip")
unzip("N03-150101_GML.zip")

ダウンロードしてきたshapefileを読み込みます。shapefileをはじめとした地理空間データをRで扱う際はsfパッケージを使うと大変楽です。

library(sf)
# Linking to GEOS 3.4.2, GDAL 2.1.2, proj.4 4.9.1

dfs.jp <- st_read("N03-20150101_GML")
# Reading layer `N03-15_37_150101' from data source `N03-20150101_37_GML/N03-15_37_150101.shp' using driver `ESRI Shapefile'
# converted into: POLYGON
# Simple feature collection with 669 features and 5 fields
# geometry type:  POLYGON
# dimension:      XY
# bbox:           xmin: 133.4465 ymin: 34.01229 xmax: 134.4475 ymax: 34.56485
# epsg (SRID):    NA
# proj4string:    +proj=longlat +ellps=GRS80 +no_defs

find_city()という関数を書きました。といってもほとんどnozmaさんのものと変わりません。対象のポリゴンと任意のポイントを指定し、ポイントがポリゴンに含まれるかを判定します。含まれる場合には市区町村コードと市区町村名を返します。

#' @param sp_polygon object of class sf, sfc or sfg
#' @param lon longitude
#' @param lat latitude
find_city <- function(sp_polygon = df, lon = lon, lat = lat) {
  
  which.row <- sf::st_contains(sp_polygon, sf::st_point(c(lon, lat)), sparse = FALSE) %>%  
    grep(TRUE, .)
  
  if (identical(which.row, integer(0)) == TRUE) {
    message("指定した座標がポリゴンに含まれません")
  } else {
  geos <- sp_polygon[which.row, ] %>%
    dplyr::mutate_if(is.factor, as.character) %>% 
    dplyr::mutate_at(dplyr::vars(N03_001:N03_004), dplyr::funs(dplyr::if_else(condition = is.na(.), true = "", false = .)))
  
    res <- tibble::data_frame(
    city_code = geos$N03_007,
    city_name = paste0(geos$N03_001, geos$N03_002, geos$N03_003, geos$N03_004)
  )
  
  return(res)
  }
}

find_city(dfs.jp, lon = 135.728965, lat = 35.039200)

## although coordinates are longitude/latitude, it is assumed that they are planar

## # A tibble: 1 × 2
##   city_code        city_name
##       <chr>            <chr>
## 1     26101 京都府京都市北区

find_city(dfs.jp, 139.884678, 35.626065)

## although coordinates are longitude/latitude, it is assumed that they are planar

## # A tibble: 1 × 2
##   city_code    city_name
##       <chr>        <chr>
## 1     12227 千葉県浦安市

find_city(dfs.jp, 139.523620, 35.495978)

## although coordinates are longitude/latitude, it is assumed that they are planar

## # A tibble: 1 × 2
##   city_code          city_name
##       <chr>              <chr>
## 1     14112 神奈川県横浜市旭区

一部の座標に関しては元のshapefileで抜け? があるためか失敗します。

find_city(dfs.jp, 139.7211411, 33.0182349)
# 指定した座標がポリゴンに含まれません
find_city(dfs.jp, 136.8231215, 34.8323039)
# 指定した座標がポリゴンに含まれません

system.time(find_city(dfs.jp, lon = 142.191237, lat = 27.0427396))
#    user  system elapsed 
#   3.676   0.223   3.917

全国のポリゴンを利用しているためか、やや時間がかかります...。市区町村名までの逆ジオコーディング手法として利用できるかとも考えましたが、件数の多い場合は別のサービス等を使うのが良さそうです。また、市区町村までいらないよ、という場合はnozmaさんの方法の通りsf::st_combine()によるポリゴンの結合を行うと良いでしょう。

jpndistrictパッケージを使う場合 (sf未対応・都道府県から検索)

座標の位置する都道府県までは判明している場合、jpndistrictパッケージを使うともう少し素早く検索できます。

jpndistrict_find_city <- function(sp_polygon = df, lon = lon, lat = lat) {
  
  which.row <- sf::st_contains(sp_polygon, sf::st_point(c(lon, lat)), sparse = FALSE) %>%  
    grep(TRUE, .)
  
  if (identical(which.row, integer(0)) == TRUE) {
    message("指定した座標がポリゴンに含まれません")
  } else {
  geos <- sp_polygon[which.row, ]
  
    res <- tibble::data_frame(
    city_code = geos$city_code,
    city_name = paste0(geos$pref_name, geos$city_name_full)
  )
  
  return(res)
  }
}

jpndistrict::spdf_jpn_pref(26) %>% 
  st_as_sf(crs = 4326) %>% 
  jpndistrict_find_city(lon = 135.728965, lat = 35.039200)

## # A tibble: 1 × 2
##   city_code         city_name
##      <fctr>             <chr>
## 1     26101 京都府京都市 北区

jpndistrict::spdf_jpn_pref(12) %>% 
  st_as_sf(crs = 4326) %>% 
  jpndistrict_find_city(139.884678, 35.626065)

## # A tibble: 1 × 2
##   city_code    city_name
##      <fctr>        <chr>
## 1     12227 千葉県浦安市

system.time(jpndistrict::spdf_jpn_pref(33) %>% 
  st_as_sf(crs = 4326) %>% 
  jpndistrict_find_city(lon = 133.93602, lat = 34.665193))
#            user          system         elapsed 
# 0.0380000000000 0.0080000000000 0.0459999999998