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

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

## パッケージ読み込み

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

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

## データ読み込み

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

```{r}
# データ読み込み
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` に格納します。

```{r}
# 以下のようなエラーが出るので、ポリゴンから無効なデータを消しておく
# 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）だけ残し、他のカラムを削除して、その後結合します。

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

## ggplotで描画する
### もっともシンプルな描画

```{r}
p<- ggplot( site ) + geom_sf()
```

遺跡ポイントデータだけを表示します。
![01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e13fa8c3-3b3f-8499-56fb-df6599ee3ef8.png)

### 市町村界を重ねる

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

```{r}
p2 <- p + geom_sf(data=town , fill = NA)
```
![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/587bb9fd-58c9-7240-1d85-c04b1722cb64.png)

### 等高線を重ねる

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

![03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/5d3441b2-6520-fc12-ede9-157c4d6bb03f.png)

### 密度図を重ねる

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

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

![05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/3ef53c54-2c5e-5a11-2cb4-016a9759f20f.png)

## 完成図

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

```{r}
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）を克服する基礎技術となるものと考えています。
![04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/ce5b919d-8970-339d-afce-0d47da982cfe.png)