北海道埋蔵文化財包蔵地データをつかって遺跡分布図をggplot2で作成する

北海道オープンデータポータルで埋蔵文化財包蔵地データが公開されました。
このデータを使って、考古学者のための遺跡分布図の作成方法を紹介します。

パッケージ読み込み

必要なパッケージはこちらです。定番パッケージに加えて、lwgeomパッケージというものをインストールします。

# パッケージ読み込み
library(tidyverse)
library(sf)
library(ggthemes)
library(viridis)
library(lwgeom)

データ読み込み

包蔵地データはポイントデータとポリゴンデータに別れています。
ポイントデータは6750件、ポリゴンデータは5498件あります。
ほかに国土数値情報の市町村界データを読み込みます。

# データ読み込み
site_po <- st_read("包蔵地_ポイント_20190913_JGD2011.shp")
site_pl <- st_read("包蔵地_ポリゴン_20190913_JGD2011.shp")
town <- st_read("N03-19_01_190101.shp")

ポリゴン重心を算出

公開されているポリゴンデータに不備があるようで、そのままではポリゴンデータをポイントデータに変換する作業がうまくいきません。
以下の方法で、不備のあるデータを除きます。
lwgeomパッケージの st_is_valid() というコマンドを使います。

st_is_valid() は不備のある行をNAで吐いてくれるので、filter() と !is.na() を組み合わせてNA以外の行（不備のない行）を抽出します。
抽出した健全なデータをパイプで st_centroid() に送って重心点を算出しています。
ポリゴン重心のポイントデータを site_pl_po に格納します。

# 以下のようなエラーが出るので、ポリゴンから無効なデータを消しておく
# CPL_geos_op("centroid", x, numeric(0), integer(0), numeric(0),  でエラー: 
#  Evaluation error: IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4.
site_pl_po <- site_pl %>% 
  filter(!is.na(st_is_valid(site_pl, reason = TRUE))) %>%
  st_centroid()

結合

ポリゴンデータ（の重心点）とポイントデータを結合します。
ポイントデータとポリゴンデータは、カラムが共通しないので、そのままでは結合できません。
遺跡名（SiteName）だけ残し、他のカラムを削除して、その後結合します。

a<-site_po %>% select(SiteName)   # ポイントデータの遺跡名だけ選択
b<-site_pl_po %>% select(SiteName)   # ポリゴン重心データの遺跡名だけ選択
site <- rbind(a,b)   # ポイントデータとポリゴン重心データを結合

ggplotで描画する

もっともシンプルな描画

p<- ggplot( site ) + geom_sf()

遺跡ポイントデータだけを表示します。

市町村界を重ねる

fill= NA を指定しないと遺跡点が上塗りされてしまいます。

p2 <- p + geom_sf(data=town , fill = NA)

等高線を重ねる

p3 <- p2 + geom_density_2d(
    aes( x = st_coordinates(site)[,1] , y = st_coordinates(site)[,2]) ,
  ) 

密度図を重ねる

密度図の場合、geom_density2d() を指定するのではなく、stat_density2d() を指定します。
ggplotチートシートにも載っている方法です。

p4 <- p2 + 
  stat_density2d( 
      aes(x=  st_coordinates(site)[,1] , y = st_coordinates(site)[,2],
      fill=..density..),
      geom = "raster",contour = FALSE
    )

完成図

説明すべきことは多いのですが、density_2d の引数には帯域幅を調整する h というオプションがあり、ここでは 0.5 を指定しています（緯度経度系なので0.5度）。
また、R本体の kde2d 関数のオプションである n というカーネル密度の解像度を設定するオプションも使えるので、ここでは片側 200pixel となるように n=200 と指定しています。

p <- ggplot(site) +
  geom_sf() +
    stat_density2d( 
      aes(x=  st_coordinates(site)[,1] , y = st_coordinates(site)[,2],
      fill=..density..),
      geom = "raster",contour = FALSE ,n = 200 , h = c(0.5 , 0.5)
    ) +
  geom_sf(data=town , fill = NA ,size=0.2 ,colour= "grey50") +
  geom_density_2d(
    aes( x = st_coordinates(site)[,1] , y = st_coordinates(site)[,2]) ,
    size = 0.3 ,colour = "White", h = c(0.5 , 0.5)
  ) +
  geom_sf( size=0.2 ,colour = "White") +
  xlim( 139,146 ) +
  scale_fill_viridis() +
  theme_minimal() +
  theme(
    axis.title= element_blank()
  )

このような分布図から文化圏を抽出することが、考古学の分布論の基本となります。
分布が何を意味するのかは、V.G．Childの古典的な研究以来、考古学の大きなテーマです。
GISを利用した空間解析が、「空白のこわさ」（佐原真1985）を克服する基礎技術となるものと考えています。