🌏Rで地域メッシュを使った地図を描きたい
この記事はFOSS4G Advent Calendar 2015 の日目です。GISとかFOSS4Gについての知識が限りなく0に近いのですが、普段使っているRでGIS的なことをしてみたい、という話です。具体的には、地域メッシュコードに基づく地図をRで描いてみる、という話になります。
地域メッシュコードでの地図が描けると、生物の分布図や広範囲での統計データのマッピングの際に便利です。実際に
匿名知的集団ホクソエムの誓いとして、必要なものがなければ作る、という言葉がありますので、それに従ってメッシュコードのデータを作成したいと思います。地域メッシュコードには、1次メッシュから順に約1km2刻みの3次メッシュまでがありますが、順にやっていきます。間に合わなかったので2次メッシュまででご勘弁ください。後でこっそり追加します
書いたRのコードは汚いので恥ずかしいため、最後に載せます。結果からご覧ください。
leafleatを使ってグリグリ動かす地図 ->
国土をカバーするメッシュ(枠)
ただの枠です。このままではよくわかりませんね。
1次メッシュ
先ほどの枠から、1次メッシュの範囲となっている部分のみを選んでいきましょう。 http://www.gsi.go.jp/MAP/HISTORY/5-25-index5-25.html このページに20万分1地勢図の区画名(及び1次メッシュコードを含んだリンク)がありますので、そこから区画名と1次メッシュコードをスクレイピングします。
先ほどの枠と組み合わせて、20万分1地勢図に記載があるメッシュのみを抽出すれば、1次メッシュコードの地図が完成です。
2次メッシュ
黒くなっているのはメッシュが集中しているせいなのですが、ちょっとずれています...(おい)。後で直します...。
3次メッシュ
計算が終わりませんでした。
R コード
library(celestial) library(broom) library(rvest) library(rgdal) library(readr) library(geojsonio) library(ggplot2) library(pforeach) # devtools::install_github("hoxo-m/pforeach") library(dplyr)
# 緯度経度の範囲 lon_range <- c(122, 154) lat_range <- c(20, 46) # 1次メッシュの間隔 (int1_lat <- dms2deg("+0d40m00.00s", sep = "dms")) # latitude
## [1] 0.6666667
(int1_lon <- dms2deg("+1d00m00.00s", sep = "dms")) # longitude
## [1] 1
seq1_lon <- seq(lon_range[1], lon_range[2], int1_lon) seq1_lat <- seq(lat_range[1], lat_range[2], int1_lat) # # 2次メッシュの間隔 (int2_lat <- dms2deg("+0d05m00.00s", sep = "dms")) # latitude
## [1] 0.08333333
(int2_lon <- dms2deg("+0d07m30.00s", sep = "dms")) # longitude
## [1] 0.125
seq2_lon <- seq(lon_range[1], lon_range[2] + int2_lon, int2_lon) seq2_lat <- seq(lat_range[1], lat_range[2] + int2_lat, int2_lat)
pforeach(i = seq1_lon, .c = rbind)( { data_frame(long = c(i, i, i + int1_lon, i + int1_lon, i)) }) -> df_long df_tmp <- data_frame(lat = rep(c(seq1_lat[2], seq1_lat[1], seq1_lat[1], seq1_lat[2], seq1_lat[2]), length(seq1_lon))) %>% bind_cols(df_long) df_res <- pforeach(i = 1:39, .c = rbind)({ dplyr::mutate(df_tmp, lat = lat + int1_lat * i) %>% rbind(df_tmp) }) %>% dplyr::mutate(group = rep(1:(nrow(.) / 5), each = 5)) df_tmp <- df_res %>% dplyr::group_by(group) %>% dplyr::summarise(code2 = as.numeric(min(lat) * 1.5), code4 = min(long) %>% substr(2, 3) %>% as.numeric(), code_1_20 = paste0(code2, code4)) df_res %<>% inner_join(df_tmp) %>% dplyr::select(long, lat, code_1_20) %>% dplyr::rename(group = code_1_20) %>% dplyr::filter(group <= 6854)
names_1_20 <- vector() code_1_20 <- vector() url <- "http://www.gsi.go.jp/MAP/HISTORY/5-25-index5-25.html" %>% read_html() html_nodes(url, css = "map area") %>% { names_1_20 <<- html_attr(., "alt") %>% Nippon::kana2roma(.) code_1_20 <<- html_attr(., "href") %>% sub(pattern = "/MAP/HISTORY/5-25-", replacement = "", x = .) %>% sub(".html$", "", .) } df_mesh <- data_frame(names_1_20, code_1_20) rm(list = grep("20$|url", ls(), value = TRUE)) dplyr::filter(df_res, group %in% df_mesh$code_1_20) %>% geojson_json(input = ., group = "group", lat = "lat", lon = "long", geometry = "polygon") %>% geojson_write(input = ., file = "jpmesh_1.geojson", # group = "group", geometry = "polygon")
リンク
過去の遺産など