読者です 読者をやめる 読者になる 読者になる

まだ厨二病

RとPyhtonを使ったデータ分析・統計解析らへんの話題をしていくだけ

🌏Rで地域メッシュを使った地図を描きたい

rgdal geojsonio celestial WIP ggplot2 mapping broom

この記事はFOSS4G Advent Calendar 2015 の日目です。GISとかFOSS4Gについての知識が限りなく0に近いのですが、普段使っているRでGIS的なことをしてみたい、という話です。具体的には、地域メッシュコードに基づく地図をRで描いてみる、という話になります。

地域メッシュコードでの地図が描けると、生物の分布図や広範囲での統計データのマッピングの際に便利です。実際に

匿名知的集団ホクソエムの誓いとして、必要なものがなければ作る、という言葉がありますので、それに従ってメッシュコードのデータを作成したいと思います。地域メッシュコードには、1次メッシュから順に約1km2刻みの3次メッシュまでがありますが、順にやっていきます。間に合わなかったので2次メッシュまででご勘弁ください。後でこっそり追加します

書いたRのコードは汚いので恥ずかしいため、最後に載せます。結果からご覧ください。

leafleatを使ってグリグリ動かす地図 ->

RPubs - メッシュ地図のテスト

国土をカバーするメッシュ(枠)

ただの枠です。このままではよくわかりませんね。

f:id:u_ribo:20151210054824p:plain

1次メッシュ

先ほどの枠から、1次メッシュの範囲となっている部分のみを選んでいきましょう。 http://www.gsi.go.jp/MAP/HISTORY/5-25-index5-25.html このページに20万分1地勢図の区画名(及び1次メッシュコードを含んだリンク)がありますので、そこから区画名と1次メッシュコードをスクレイピングします。

先ほどの枠と組み合わせて、20万分1地勢図に記載があるメッシュのみを抽出すれば、1次メッシュコードの地図が完成です。

f:id:u_ribo:20151210055206p:plain

2次メッシュ

黒くなっているのはメッシュが集中しているせいなのですが、ちょっとずれています...(おい)。後で直します...。

f:id:u_ribo:20151210055220p:plain

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")

リンク

過去の遺産など

広告を非表示にする