recogeoパッケージで市町村合併に伴う再集計をやってみる
市区町村単位で集計された統計データで時点間の比較をしようとすると、市町村合併によって2時点間で地域の対応が取れないことがあります。必要に応じて自前でスクリプトを書いて対応していますが、よさげなパッケージ（recogeo package）がリリースされていたのでちょっと試してみました。
2005年と2010年の国勢調査で得られた福岡県の人口データを例にやってみます。データは統計GISの小地域境界データのShapeファイルを使います。
パッケージのロード
library(tidyverse) library(sf) library(rmapshaper)
2005年データのダウンロードと展開
url2005 <- "https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212005&code=40&coordSys=1&format=shape&downloadType=5" tf <- tempfile() curl::curl_download(url2005, tf) unzip(tf, exdir = "pop2005", junkpaths = TRUE) unlink(tf) rm(tf)
2005年データの読み込み
st_read()でシェープファイルを読み込む。海域のポリゴンを除去するために、HCODE == 8101で町丁・字等に限定。rmapshaper::ms_dissolve()で市区町村単位でディゾルブ。市区町村コードをPREFとCITYを結合して作成。rmapshaper::ms_simplify()でジオメトリを簡素化。rmapshaper::ms_filter_islands()で面積の小さい離島を除去。
pop2005 <- st_read("./pop2005", options = "ENCODING=CP932", stringsAsFactors = FALSE) %>% filter(HCODE == 8101) %>% select(PREF, CITY, CITY_NAME, JINKO, SETAI) %>% ms_dissolve("CITY_NAME", sum_fields = c("JINKO", "SETAI"), copy_fields = c("PREF", "CITY")) %>% mutate(CITY_CODE = str_c(PREF, CITY)) %>% select(-PREF, -CITY) %>% ms_simplify() %>% ms_filter_islands(min_area = 2000000) %>% rename(JINKO2005 = JINKO, SETAI2005 = SETAI) %>% st_set_crs(pop2005, 4612) %>% st_transform(2444)
2010年も同様に
url2010 <- "https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212010&code=40&coordSys=1&format=shape&downloadType=5" tf <- tempfile() curl::curl_download(url2010, tf) unzip(tf, exdir = "pop2010", junkpaths = TRUE) unlink(tf) rm(tf) pop2010 <- st_read("./pop2010", options = "ENCODING=CP932", stringsAsFactors = FALSE) %>% filter(HCODE == 8101) %>% select(PREF, CITY, CITY_NAME, JINKO, SETAI) %>% ms_dissolve("CITY_NAME", sum_fields = c("JINKO", "SETAI"), copy_fields = c("PREF", "CITY")) %>% mutate(CITY_CODE = str_c(PREF, CITY)) %>% select(-PREF, -CITY) %>% ms_simplify() %>% ms_filter_islands(min_area = 2000000) %>% rename(JINKO2010 = JINKO, SETAI2010 = SETAI) %>% st_set_crs(pop2005, 4612) %>% st_transform(2444)
地図の描画
市区町村数は市町村合併により、97から72に減少しています。
nrow(pop2005) ## [1] 97 nrow(pop2010) ## [1] 72
geom_sf()で地図の描画。
gridExtra::grid.arrange( ggplot(pop2005) + geom_sf(aes(fill = log(JINKO2005)), lwd = 0.2) + scale_fill_gradient(low="white", high="red"), ggplot(pop2010) + geom_sf(aes(fill = log(JINKO2010)), lwd = 0.2) + scale_fill_gradient(low="white", high="red"), ncol = 2)
recogeoパッケージを使う
ここからやっと本題。recogeo::reconcileGeographies()関数を使うと、2つの空間オブジェクトに含まれるポリゴン同士の関係（同じsame / 含むcontains / 重なるintersects）の判定をしてくれます。ただ、実際にはこれらの判定がそのままではうまくいかないことがあります。例えば、若松区（40103）の場合、
pop2005.wakamatsu <- pop2005 %>% filter(CITY_NAME == "若松区") pop2010.wakamatsu <- pop2010 %>% filter(CITY_NAME == "若松区") ggplot() + geom_sf(data = pop2005.wakamatsu, color = "red", fill = NA) + geom_sf(data = pop2010.wakamatsu, color = "blue", fill = NA)
このように微妙にずれていますので、実際には同じ自治体なのですが、いずれも包含関係になりません。
そこで片方について、500mバッファをとってみます。
ggplot() + geom_sf(data = pop2005.wakamatsu, color = "red", fill = NA) + geom_sf(data = st_buffer(pop2010.wakamatsu, 500), color = "blue", fill = NA)
これなら確実に2010年のポリゴンで2005年のポリゴンを含むことができます。逆も同様に実行した場合に、包含関係が成立すれば、2つのポリゴンは同一自治体のものと判定するという感じになっている（と推測してます）。一方の包含関係のみが成立する場合はcontainになって、今回のデータの場合は、合併があったことを示すものになります。
次に、intersectについてですが、そのままintersectかどうかのチェックをすると、多くの隣接自治体がintersectと判定されてしまいます。2005年の若松区は2010年の若松区と芦屋町と八幡西区とintersectの関係にあると判定されます。
st_intersects(pop2005.wakamatsu, pop2010) ## 1: 3, 8, 49 ggplot() + geom_sf(data = pop2010[c(3, 8, 49),], color = "blue", fill = NA) + geom_sf(data = pop2010.wakamatsu, color = "red", fill = NA)
これを抑止するために、intersectしたときの共通部分の面積が、一定以上の値でないと、intersectしているとみなさないということにします。今回のデータでは合併のみしか起こっていませんので、intersectの判定は生じないはずです。つまり、intersectが出現しないように、その閾値を設定すればよいことになります。
実際の判定に、recogeo::reconcileGeographies()関数を使います。
devtools::install_github('fraba/recogeo') library(recogeo) pop_recogeo <- reconcileGeographies( pop2010, pop2005, "CITY_CODE", "CITY_CODE", dist_buffer = 500, min_inters_area = 30000)
これで、pop_recogeoにこの期間での合併情報がdata.frameの形式で格納されます。以下は、一部の出力ですが、2006年にあった飯塚市（40205）と周辺の4自治体の合併情報が確認できます。
pop_recogeo %>% arrange(unigeokey_A) %>% slice(18:24) ## unigeokey_A unigeokey_B type ## 1 40205 40205 AcontainsB ## 2 40205 40425 AcontainsB ## 3 40205 40426 AcontainsB ## 4 40205 40427 AcontainsB ## 5 40205 40428 AcontainsB ## 6 40206 40206 same ## 7 40207 40207 same
この情報を使って、recogeo::reconcileData()関数で人口と世帯数の集計をしてみます。
pop_recodat <- reconcileData(pop_recogeo, pop2010, pop2005, "CITY_CODE", "CITY_CODE", varA = c("JINKO2010", "SETAI2010"), varB = c("JINKO2005", "SETAI2005"), return_spatial = "A") pop_recodat %>% as_tibble() %>% select(1:5) %>% slice(1:5) ## # A tibble: 5 x 5 ## .unigeokey_new JINKO2010 SETAI2010 JINKO2005 SETAI2005 ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 1 104469 44232 108677 44364 ## 2 2 85167 33495 87340 32754 ## 3 3 61583 28064 63714 28417 ## 4 4 181936 89036 183286 87459 ## 5 5 214793 86420 214624 83716
2010年の行政区域で2005年の人口と世帯数が再集計されています。 これで、世帯当たりの人数を地図上に可視化してみると次のようになります。
pop_recodat %>% mutate(`2010` = JINKO2010 / SETAI2010, `2005` = JINKO2005 / SETAI2005) %>% select(`2005`, `2010`) %>% gather(year, n_household, `2005`, `2010`) %>% ggplot(aes(fill = n_household)) + geom_sf() + facet_wrap(.~year) + scale_fill_gradient(low="white", high="red")
とりあえず、合併のみという単純な場合なら使えそうな感じです。政令指定都市への移行などで生じる分割や、もっと複雑な場合もありますが、その検証はまた別の機会に。
Enjoy!
