sfパッケージ&国土数値情報でコロプレス図
Rで地図を描画する際はspパッケージが必須だったわけですが、今後sfパッケージの方が便利にになるらしいです。dplyrと相性がいいんでしょうね。この本でやっている、ごみ排出量のコロプレス図の作成をsfパッケージでやってみることにします。
シェープファイルの取得
国土数値情報からシェープファイルをダウンロード。kokudocuuchiパッケージを使って、URLを確認。
library(dplyr) kokudosuuchi::getKSJSummary() %>% filter(title=="行政区域")
## # A tibble: 1 x 5 ## identifier title field1 field2 areaType ## <chr> <chr> <chr> <chr> <chr> ## 1 N03 行政区域 政策区域 - 3
ということなので、福岡県(prefCode = 40)を指定してURL一覧を取得。そこから2015年のもののみ取り出す。
kokudosuuchi::getKSJURL("N03", prefCode = 40) %>% filter(year==2015) %>% select(zipFileUrl) %>% pull(1) -> url
ファイルをダウンロードしてカレントディレクトリに展開。
tf <- tempfile() curl::curl_download(url, tf) unzip(tf, exdir = getwd(), junkpaths = TRUE) unlink(tf) rm(tf)
シェープファイルの読み込み→加工→白地図の描画
sfパッケージをロードして、シェープファイルを読み込む。
library(sf) d <- st_read(".", options = "ENCODING=CP932", stringsAsFactors = FALSE)
政令指定都市については、ごみ排出量のデータは、全体での集計になっているのに対して、シェープファイルは区ごとのデータになっているため、政令指定都市全体をaggregate()関数を使って1つのポリゴンにまとめる。まとめた後に、市区町村名(N03_004)と、市区町村コード(N03_007)の修正が必要。
d %>% filter(N03_003 %in% c("北九州市","福岡市")) %>% aggregate(by = list(.$N03_003), FUN = "head", n=1) %>% select(-1) %>% mutate(N03_004 = N03_003, N03_007 = case_when(N03_004 == "北九州市" ~ "40100", N03_004 == "福岡市" ~ "40130")) -> d.seirei
他の市区町村についても、国土数値情報のデータは1つの市町村が多くのポリゴンの集まりになっているので、これらも1つのポリゴンにまとめる。
d %>% filter(!(N03_003 %in% c("北九州市", "福岡市"))) %>% aggregate(by = list(.$N03_007), FUN = "head", n=1) %>% select(-1) -> d.others
これらをまとめる。
rbind(d.seirei, d.others) %>% arrange(N03_007) -> d.new
白地図を出力。
plot(st_geometry(d.new), border="grey", col="white")
ごみ排出量データの取得と、sfオブジェクトへのマージ
ダウンロードするファイルのURLと保存先ファイル名の指定。e-statのAPI使ってやると便利かもしれないけど、idの取得が必要になるので今回はパス。
waste_url <- "http://www.env.go.jp/recycle/waste_tech/ippan/h27/data/shori/city/40/01.xls" waste_fn <- "waste_fuk_h27.xls"
ダウンロードして、readxlパッケージでExcelのファイルを読み込む。列名を設定して、特に使う項目は分かりやすい列名に変更。
curl::curl_download(waste_url, waste_fn) waste_fuk <- readxl::read_excel(waste_fn, sheet = 1, skip = 6) colnames(waste_fuk) <- paste0("V", seq(1:ncol(waste_fuk))) waste_fuk %>% rename(jcode=V2, name=V3, waste_total=V12) -> waste_fuk
sfオブジェクトへマージ。
left_join(d.new, waste_fuk, by = c("N03_007" = "jcode")) -> d.new2
コロプレス図(塗り分け地図)の作成。
library(classInt) library(RColorBrewer) ci <- classIntervals(d.new2$waste_total, style = "pretty") d.new2 %>% select(waste_total) %>% plot(col=findColours(ci, rev(brewer.pal(7, "RdYlGn"))))
ggplot2のgeom_sf()を使った地図の描画
開発版のggplot2では、geom_sf()を使って、sfオブジェクトからggplot2形式のプロットが作成できます。 Windowsの場合は開発ツールをインストールしてから、以下のようにすればインストールできます。
devtools::dev_mode() devtools::install_github("hadley/scales") devtools::install_github("tidyverse/ggplot2")
library(ggplot2) ggplot(d.new2, aes(fill=waste_total)) + geom_sf() + scale_fill_distiller(palette="RdYlGn") + theme_bw()
enjoy!