Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

@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 
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away