홋카이도 매장 문화재 포장지 데이터를 이용해 유적 분포도를 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)을 극복하는 기초 기술이 될 것으로 생각하고 있습니다.