Last active
April 30, 2017 10:40
-
-
Save nozma/3c564a2cc8eec9ac340954df766762c9 to your computer and use it in GitHub Desktop.
Rで緯度経度から都道府県を求める ref: http://qiita.com/nozma/items/808bce2f496eabd50ff1
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| library(sf) | |
| # 日本地図のシェープファイルを読み込む | |
| f <- "N03-160101_GML/N03-16_160101.shp" | |
| japan <- st_read(f) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| > 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... |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # 都道府県名一覧の取得 | |
| 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より早い |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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) | |
| } |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| > 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 |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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))),] | |
| } |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| > 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... |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| > 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 |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #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))),] | |
| } |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| > system.time(find_pref(p1)) | |
| ユーザ システム 経過 | |
| 0.003 0.001 0.004 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment