cucumber flesh

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

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

追記 (180814)

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!

geofacetで日本のデータが利用可能になりました

もう数ヶ月も前の話ですが、CRANに登録されているgeofacetパッケージが更新されて、バージョンが0.1.9となりました。geofacetパッケージは、ggplot2ベースの作図パッケージの一種ですが、グラフを地理的な空間関係と対応付けて描画するという特徴があります。例えば、次のような図を見たことがある人がいるのではないでしょうか。

f:id:u_ribo:20180622064350p:plain

この図はアメリカでの2000年から2016年にかけての失業率を州に分けて表示しています。ただ図を並べるだけでなく、各州の空間関係を考慮した配置をすることで各州の関係がわかりやすく可視化されていると思います。この図はgeofacetのサンプルコードとして掲載されているもので、下記を実行することで描画されます。

library(ggplot2)
library(geofacet)

ggplot(state_unemp, aes(year, rate)) +
  geom_line() +
  facet_geo(~ state, grid = "us_state_grid2") +
  scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  ylab("Unemployment Rate (%)")

geofacetには、こうした図を書くために、いくつかの国や地域に対応していますが、バージョン0.1.9で日本にも対応するようになりました!追加した当人として、geofacetで日本のデータを利用する例を紹介します。アメリカに比べると少し不恰好で、現実の大きさや配置を再現できていない部分もあるかと思いますが、ぜひお使いください!

github.com

続きを読む

ggplot2のsizeが意味するもの

付け焼き刃な知識故に、間違いございましたら指摘頂けますと幸いです。

ggplot2を使っていると、数種類の“size”を設定するタイミングがある。例えば、geom_point()で散布図を描画する際のポイントの大きさを指定するためのsize、または連続した線を描画するのに用いられるgeom_line()でのsize、あるいは文字を重ねるためのgeom_text()annotate()の中でのsizeである。また、少し異なる引数名であるがテーマ指定の関数theme_*()にも“base_size”という名の引数が用意されている。ここで肝心なことは、テキストとその他では、同じsizeであっても単位が異なるということだ。このことは次のStackOverflowの投稿にまとめられている。

stackoverflow.com

まとめて書いてしまうとこの通り。

  • レイヤーとして重ねるgeom_point()geom_line()などの要素のサイズはmm単位(ポイントの大きさ、線の太さなど)
  • 軸やタイトルなど、フォントを表示するテキスト系のサイズはpoint単位(1pointsはおよそ0.352778mm)
  • geom_point()のポイントの輪郭の太さを制御するstrokeの単位もpoint(ggplot2::.stroke

以下、この記事ではそのことを確かめてみたいと思う。

適当なggオブジェクトを用意し、それに手を加えて比較していくことにする。

library(ggplot2)
# 日本語の表示のためにフォントを指定
theme_set(theme_gray(base_family = "IPAexGothic"))

p <- 
  ggplot(mtcars, aes(wt, mpg))

まずはgeom_point()geom_line()でのsizeである。ここではsize = 10を指定した(既定では1.5)。

p_out <- 
  p + 
  geom_point(size = 10)

この結果をggsave()で出力した図を表示する。ポイントのサイズが10mmにしたので、検証のために出力時のサイズはその10倍となる10cmとした。

ggsave(filename = "out.png", 
       plot = p_out, 
       width = 10, 
       height = 10, 
       units = "cm")

f:id:u_ribo:20180611230853p:plain

確かにmmらしい。もっとも、これはソースコードにも示されている。

Size examples Should be specified with a numerical value (in millimetres), or from a variable source

github.com

次に、軸やタイトルテキストのsizeである。軸のラベルを除いて、size = 10を指定した。

p2 <- 
  p_out +
  # 適当な場所に文字「A」を表示
  annotate("text", x = 2, y = 15, size = 10, label = "A", color = "red") +
  ggtitle("文字の大きさはpointで指定")

p2 + 
  theme(title = element_text(size = 10),
         axis.title = element_text(size = 10),
         # グラフの両軸に用いられる数値の大きさ
         axis.text = element_text(size = 5))

f:id:u_ribo:20180611231008p:plain

フォントを扱っていて、なおかつ同じsize = 10を指定していてもannotate()が表示するフォントが大きい。これはannotate()がレイヤーの要素として評価されるためで、geom_point()のようにsizeはmm単位での指定となるためだ。

もし、タイトルや軸のサイズをmm単位にしたい時はpointをmmへ変換する必要がある。ggplot2::.ptはggplot2の内部で用いられている変数で、ドキュメントもある

p2 +
  theme(axis.title = element_text(size = ggplot2::.pt * 10))

f:id:u_ribo:20180611231028p:plain

最後にtheme_*()の引数“base_size”だが、こちらも単位はpointである。すなわち次の2つのコードによって出力される図は等しい。

ggplot2が生成するような描画オブジェクトの差分を確認するにはvdiffrが便利である。ここではwidget_diff()を用いて差分の検証を行うが、GitHubのように2つの画像を表示し、スライドさせながら差分を確認するwidget_slide()やテスト時に役立つ関数も備わっている。

qiita.com

# plotとaxisのラベルの大きさは、テーマのbase_sizeに対してI()で指定した値をかけた大きさとなる
p_default <- 
  p2 + 
  theme_gray(base_size = 11,
             base_family = "IPAexGothic")

p_custom <- 
  p2 +
  theme_replace(
    plot.title = element_text(size = 11 * I(1.2)),
    axis.title = element_text(size = 11 * I(0.8)))
# 2つのオブジェクトで差分がない場合は白い画面が表示される
vdiffr::widget_diff(
  p_default,
  p_custom
)

# 変更してみる
# フォントサイズが変更されたことで全体の配置も微妙にずれる
vdiffr::widget_diff(
  p_default,
  p_custom +
    theme(plot.title = element_text(size = 20))
)

f:id:u_ribo:20180611231716p:plain

学術誌に投稿する図では、軸のテキストサイズやらについて細かく指定しされているものもある。怒られないように、少しずつggplot2に慣れていこうと思う。 Enjoy!