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

北海道オープンデータポータル 埋蔵文化財包蔵地データをつかったカーネル密度推定の手法を紹介します。

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

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

# パッケージ読み込み
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 に格納します。

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

結合

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

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

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

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

カーネル密度推定の実行

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

# 遺跡一データを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をコピペしました。

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をコピーするのが便利です。

出力したgeotiffをQGISで表示

ggplotで描画したい

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

# データフレームに変換
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")  # カラム名を付値

# 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()
   )

まとめ

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