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

f:id:nonki1974:20170712081741j:plain

ごみ排出量データの取得と、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"))))

f:id:nonki1974:20170712082011j:plain

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

f:id:nonki1974:20170712082145j:plain

enjoy!