cucumber flesh

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

st_union()で埋まらないポリゴンの穴を埋める

sf sfheaders

追記 (2020-03-21)

sfheadersパッケージ (v0.2.1) で実装された sf_remove_holes() を使っても隙間を埋められます。

library(sfheaders)
sf_pref33_fill <- 
  sf_pref33 %>% 
  st_union() %>% 
  sf_remove_holes()

追記 (2018-08-14)

id:yutannihilation さんから、Twitter経由でrmapshaperパッケージのms_dissolve関数を使う方法を教えてもらいました。この方法も簡単です。どうもありがとうございます！

はじめに

この記事の中では、国土数値情報が提供する「行政区域データ」を例にしていますが、最新のデータとなる「平成30(2018)年1月1日時点」のものを使えばこのような問題は生じません。もし過去の行政区域データや同様の問題に直面した際にご参考ください。

問題

複数のポリゴンをsf::st_union()を使って一つのポリゴンに結合しようとするとき、ポリゴンの中にわずかな溝のようなものが生成されてしまうことがあります。こんな感じで。

f:id:u_ribo:20180813234741p:plain

これは対象にしているポリゴンが接する面の間にわずかに隙間があるために生じる問題です。今日はこうしたポリゴンの隙間を綺麗に埋める方法を紹介します。まずはデータを用意しましょう。

行政区域データの用意

国土数値情報の「行政区域」を利用します。ここでは岡山県のデータを用意することにしました。以下のコードでデータのダウンロードもRから実行します。ダウンロード先のURLはkokudosuuchiを使って求めます。岡山県の都道府県コードは33であり、データ抽出の処理として実行するdplyr::filter()の中でareaCodeとしてその値を与えます。

library(dplyr)
library(sf)

# 行政区域データのデータ一覧を取得
df_source <-
  kokudosuuchi::getKSJURL("N03")

# ダウンロード先のURL
urls <-
  df_source %>%
  filter(year == 2017, areaCode == 33) %>%
  magrittr::use_series(zipFileUrl)

行政区域のデータは、いくつかの時期に渡って提供されていますが、データの基準となっているのが2017年のものを用意しました。繰り返しとなりますが、ここでは例のために最新でない2017年のデータを使います。最新の2018年データでは、岡山県に関わらず全国で今回のような問題が発生しません。

# 圧縮ファイルをダウンロード。解凍して圧縮ファイルは削除する
download.file(urls, basename(urls))

files <-
  list.files(pattern = "N03-170101_33", full.names = TRUE)

unzip(files, 
      stringr::str_remove(basename(files), ".zip"))
unlink(files)

問題の再現: st_union()で埋まらない隙間

ダウンロードしたshapeファイルをR上に読み込んで、ポリゴンをプロットしてみます。

sf_pref33 <- 
  st_read(here::here("N03-170101_33_GML"), 
          stringsAsFactors = FALSE, 
          quiet = TRUE,
          options = c("ENCODING=CP932"))

綺麗な岡山県のポリゴンですね。

plot(st_geometry(sf_pref33))

f:id:u_ribo:20180813234847p:plain

sfオブジェクトを可視化する方法については以前の記事も参考ください。

uribo.hatenablog.com

結合する

それでは、問題の結合処理を行いましょう。ポリゴンを結合するのにst_union()を使うかと思いますが、ここでst_union()との違いを整理しておきます。2つの関数は似ていますが、下記の点で異なります。

  • st_combine()... 境界を維持したまま一つのgeometryにまとめる（複数のPOLYGONがMULTIPOLYGONになる）。
  • st_union()... 与えられたgeometryを一つのgeometryにする。geometry同士が重なる境界がある場合にその境界を削除する。

今回は、市区町村単位での行政区域に分かれているポリゴンデータを結合して、一つの「岡山県」ポリゴンにしたいのでst_union()を使うことになります。そして関数を適用すると冒頭の穴あき問題に衝突します。

sf_pref33_union <- 
  sf_pref33 %>% 
  # by_feature = TRUEにすると、データフレームの列 (feature)の要素ごとの結合になるのでFALSEを指定します（既定値）
  st_union(by_feature = FALSE) %>% 
  st_transform(crs = 4326)

結合したポリゴンを確認して見ると、わずかに空隙があることがわかります。st_union()の境界線の結合は、ポリゴンが接する際に行われるので、元のポリゴンにわずかでも隙間があると、その隙間が維持されてしまいます。

# 冒頭の図を再現するコード
# 図の出力は省略します
library(ggplot2)

ggplot() +
  geom_sf(data = sf_pref33_union, fill = "transparent")

なお、st_uinon(by_feature = FALSE)を実行した時の返り値はsfcなのでsfを維持したい時はsummarise()を使うと良いですが、ポリゴンの穴問題は解決されるわけではありません。

ggplot() +
  geom_sf(data = sf_pref33 %>% 
            summarise(prefecture = "岡山県"), 
          fill = "transparent")

解決策

といっておきながら、GitHub Issuesに解決策がまとめられています。その方法を真似してみましょう。

github.com

sf_pref33_union_fix <- 
  # st_union()を実行し、MULTIPOLYGONとなったPOLYGONを元に戻します。
  # featureの情報は元に戻らないので注意です
  sf_pref33_union %>% 
  st_cast("POLYGON") %>% 
  purrr::map(
    ~ .x[1]
  ) %>% 
  # ポリゴンを形成する座標からMULTIPOLYGONを作成
  st_multipolygon() %>% 
  # sfcに変換
  st_sfc(crs = 4326)

確認します。おお、できました。これが求めていたものです。

ggplot() +
  geom_sf(data = sf_pref33_union_fix, fill = "transparent")

f:id:u_ribo:20180813235059p:plain

視覚による確認に加えて、失敗と成功の2つのgeometryのうち、重複のない部分をst_difference()で抽出してみましょう。この処理の中で穴あきポリゴンの方にlwgeom::st_make_valid()を使っていますが、これはgeometryがinvalidであるために怒られるので、それを回避するために、つまりvalidなポリゴンに修正するためです（この時点でちょっと変なポリゴンだということがわかる）。

# 時間がかかります
st_pref33_diff <- 
  st_difference(sf_pref33_union_fix, 
                lwgeom::st_make_valid(sf_pref33_union))

ggplot() +
  geom_sf(data = st_pref33_diff)

f:id:u_ribo:20180813235137p:plain

問題の部分が取り出せたと思います。

最初私は、各geometryに対してst_buffer()を適用し、隙間を埋めれば良いのでは？と思いましたが、それだと外側が膨らんでしまうのでだめでした。また、特定のgeometryだけに適用するのも、どのgeometryが原因となっているのかを探る必要があったりとうまくいきませんでした。というわけで最初の方法を使うのが良いかと思います。

他に良い方法があれば教えてください。Enjoy!