고고학을 위한 커널 밀도 추정

홋카이도 오픈 데이터 포털 매장 문화재 포장지 데이터를 사용한 커널 밀도 추정의 수법을 소개합니다.

홋카이도 매장 문화재 포장지 데이터를 이용해 유적 분포도를 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 클래스로 출력되기 (위해)때문에 취급하기 쉽다고 생각합니다.