考古学のためのカーネル密度推定

[北海道オープンデータポータル](https://www.harp.lg.jp/opendata/) 埋蔵文化財包蔵地データをつかったカーネル密度推定の手法を紹介します。

[北海道埋蔵文化財包蔵地データをつかって遺跡分布図をggplot2で作成する](https://qiita.com/ishiijunpei/items/7e934eb1e96f0fe0bf50) ではggplot2の機能を利用したカーネル密度推定図を作成しましたが、より汎用性が高く、使い回しの効くカーネル密度推定図を作成します。

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

```{r}
# パッケージ読み込み
library(tidyverse)
library(sf)
library(ggthemes)
library(viridis)
library(lwgeom)
library(GISTools)
library(raster)
# データ読み込み
# 包蔵地ポイントデータ
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` は不備のあるデータをNAで吐いてくれるので、filter() と !is.na() を組み合わせてNA以外の行（不備のない行）を抽出します。
抽出した健全なデータをパイプで `sf::st_centroid()` に送って重心点を算出しています。
ポリゴン重心のポイントデータを `site_pl_po` に格納します。

```{r}
site_pl_po <- site_pl %>% 
  filter(!is.na(st_is_valid(site_pl, reason = TRUE))) %>% # NA行以外の行を検索
  sf::st_centroid() # ポリゴン重心を算出
```

## 結合

ポリゴンデータ（の重心点）とポイントデータを結合します。

```{r}
# 遺跡名（SiteName）のみを抽出
a<-site_po %>% dplyr::select(SiteName)  # ポイントデータ
b<-site_pl_po %>% dplyr::select(SiteName) # ポリゴン重心データ
site <- rbind(a,b)
```

ついでに、市町村界も北海道全体を一つのポリゴンにまとめてしまいます。

```{r}
# 市町村を結合
hkd <- sf::st_union(town)
```

## カーネル密度推定の実行

`GISTools::kde.points()`でカーネル密度を算出します。
出力は`SpatialPixelsDataFrame`なので、geotiffに書き出したり、ラスタ演算に利用できます。
`kde.points()`で指定する引数 n=グリッドの数（解像度）、h = カーネル密度推定の半径 です。

```{r}
# 遺跡一データをspクラスに変換してkde.points
dens <- site %>% 
  as("Spatial") %>%   #spクラスに変換
  GISTools::kde.points(n=500 , h = 0.5)  # 片側500ピクセル　検索半径0.5度
```

## geotiffに出力
カーネル密度をgeotiffに出力します。
投影系で出力したかったので、`spTransform`で投影変換します。
指定する`CRS=`はQGISのプロジェクトのCRSをコピペしました。

```{r}
dens %>% 
  spTransform(
    CRS = "+proj=utm +zone=54 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") %>% 
  raster()  %>% 
  writeRaster("dens.tif" , overwirte = TRUE)
```

QGISのプロジェクトのプロパティからJGD2000 UTMzone54のCRSをコピーするのが便利です。
![prj.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e76111e1-b14d-18b2-87e1-b7efdc38231e.png)


出力したgeotiffをQGISで表示
![dens_QGIS.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/3d8d4169-3c38-4695-19a0-88b50b8ebcea.png)


## ggplotで描画したい
`GISTools::level.plot()`でも描画できますが、やはりggplotで描画したいです。
手始めにデータフレームに変換します。

```{r}
# データフレームに変換
dens_d <- dens$kde %>% as.data.frame()  # 密度データ
dens_x <- dens$Var1 %>% as.data.frame() # x座標
dens_y <- dens$Var2 %>% as.data.frame() # y座標
dens_data <- dplyr::bind_cols(dens_d,dens_x,dens_y)  # 密度、x座標、y座標を結合
colnames(dens_data) <- c("dens","x","y")  # カラム名を付値
```

```{r}
# ggplotで描画 
dens_data %>% 
  ggplot() +
    geom_raster(aes(x=x,y=y,fill=dens)) +
    geom_sf(data=hkd , fill = NA ,size=0.2 ,colour= "grey50") +
    geom_sf(data= site ,size=0.2 ,colour = "White") +
    geom_contour(
      aes(x = x , y = y , z = dens) ,
      size = 0.3 ,colour = "White",
    ) +
    xlim( 139,146 ) +
   scale_fill_viridis() +
   theme_minimal() +
   theme(
      axis.title= element_blank()
   )
```
![r_density.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/2e84a69f-a862-010f-a826-8a45a701da72.png)

## まとめ
kernel密度推定を実行できる関数はsplancsパッケージやspatstatパッケージなどにも含まれていますが、`GISTools::kde.points()`がspクラスで出力されるため扱いやすいと思います。