考古学のためのRでGIS〜遺跡の標高と海岸線からの距離を分析する〜

道南ブロック博物館施設等連絡協議会のブログに「コラムリレー（第193回）不思議な旧石器時代の遺跡分布」を投稿しました。考古学で利用頻度が高そうな分析を行ったので、手順を紹介します。内容は初心者向けです。

パッケージを読み込み

利用するパッケージはsf、terra、tidyverseです。
sfはベクタデータ、terraはラスタデータの読み込みと加工に使います。

# パッケージ読み込み
library(sf)
library(terra)
library(tidyverse)

データ読込

DEMは国土地理院基盤地図情報数値標高モデル10mメッシュ、海岸線は国土交通省国土数値情報、遺跡ポイントデータは埋蔵文化財包蔵地（GISデータ）【北海道】を利用しました。

# DEM
dem <-
  terra::rast("elevation_utm54/DEM_JGD2000utm54.tif") 
# 遺跡ポイント
site <-
  sf::st_read("site_JGD2000UTMzone54.gpkg")
# 海岸線
coast <-
  sf::st_read("CostLineHKD.gpkg") %>% 
  summarise()

埋蔵文化財包蔵地データはポイントデータとポリゴンデータの2種があるのですが、ポリゴンデータはsf::st_centroid()関数で重心位置を算出してポイントデータに変換しています。

読み込んだデータを描画します。
うまく読み込めているようです。

# 遺跡表示
plot(dem)
plot(st_geometry(site), add = T)
plot(coast, add = T)

DEMの標高値をポイントデータに付値する

遺跡立地点の標高値をDEMから取得します。ラスタの値を抽出するコマンドはsfにもterraにもあるのですが、圧倒的にterraパッケージが高速です。terra::extract()コマンドは重複する位置にある他のGISデータの値を取得します。

# 遺跡立地地点の標高値を抽出
elv <-
  terra::extract(
    x = dem,
    y = vect(site) # terraのSpatVector形式に変換
  ) %>% 
  rename(Elv = DEM_JGD2000utm54) 

抽出した標高値はこのようなdata.frameに格納されます。

> head(elv)
  ID Elv
1  1  10
2  2  16
3  3   5
4  4  21
5  5  22
6  6  25

続いて、抽出した標高値を元の遺跡ポイントに結合します。あわせてdplyr::select()関数で不要な項目を削除しています。

# 遺跡ポイントに標高値を付値
site <-
  site %>% 
  bind_cols(elv$Elv) %>%
  dplyr::rename(Elv = ...23) %>% # 標高値をrename
  dplyr::select(名称, 所在地, 種別, 時代, 立地, 標高, 出土遺物, Elv)

Elvに標高が格納されました。

名称	所在地	種別	時代	Elv	geom
坊主山遺跡	江別市対雁1・2	墳墓	続縄文、擦文、アイヌ	10	POINT (543667.1 4774414)
江別第2チャシ跡	江別市対雁125-1ほか	チャシ跡	アイヌ	16	POINT (542825.3 4774351)
美原遺跡	江別市美原河川敷地	遺物包含地	不明	5	POINT (546291.3 4774576)

旧石器出土遺跡と非旧石器の遺跡を区分する

一つの埋蔵文化財包蔵地に複数の時代の遺跡が存在する場合、「時代」のフィールドには複数の時期が入力されていることがあります。以下のコードで旧石器出土遺跡と旧石器が出土しない遺跡に区分します。この場合、旧石器出土遺跡である「paleo」属性をもつ遺跡には、旧石器以外の時代が含まれることになります。

dplyr::mutate()関数で新たに「class」フィールドを作成し、stringr::str_detect()関数で「旧石器」という単語が含まれる場合は「paleo」、含まれない場合は「non-paleo」という値を入力します。

# 旧石器出土遺跡と非旧石器の遺跡の区分
site <-
  site %>% 
  dplyr::mutate(
    class = case_when(
      str_detect(時代,"旧石器") ~ "paleo", # 旧石器遺跡はpaleo
      !str_detect(時代,"旧石器") ~ "non-paleo" # 非旧石器遺跡はnon-paleo
    )
  )

名称	所在地	種別	時代	Elv	geom	class
坊主山遺跡	江別市対雁1・2	墳墓	続縄文、擦文、アイヌ	10	POINT (543667.1 4774414)	non-paleo
江別第2チャシ跡	江別市対雁125-1ほか	チャシ跡	アイヌ	16	POINT (542825.3 4774351)	non-paleo
美原遺跡	江別市美原河川敷地	遺物包含地	不明	5	POINT (546291.3 4774576)	non-paleo

海岸線からの距離を算出

sf::st_distance()関数で遺跡ポイントから海岸線までの直線距離を計算します。

# 海岸線からの距離計算
site_near <- 
  sf::st_distance(site, coast)
> head(site_near)
Units: [m]
         [,1]
[1,] 20758.26
[2,] 20135.15
[3,] 22734.56
[4,] 20824.10
[5,] 20456.33
[6,] 20408.00

遺跡データに海岸線からの距離を結合

標高を付値したのと同じやり方で、海岸からの距離を遺跡データに付値します。

site <-
  site %>% 
  bind_cols(site_near) %>% 
  dplyr::rename(coast_dist = ...8)

名称	所在地	種別	時代	Elv	class	coast_dist	geom
坊主山遺跡	江別市対雁1・2	墳墓	続縄文、擦文、アイヌ	10	non-paleo	20758.26	POINT (543667.1 4774414)
江別第2チャシ跡	江別市対雁125-1ほか	チャシ跡	アイヌ	16	non-paleo	20135.15	POINT (542825.3 4774351)
美原遺跡	江別市美原河川敷地	遺物包含地	不明	5	non-paleo	22734.56	POINT (546291.3 4774576)

海岸からの距離をヒストグラムで可視化

旧石器出土遺跡と旧石器が出土しない遺跡の海岸からの距離を比較します。

site %>% 
  drop_na(paleo) %>% 
  ggplot() +
  geom_histogram(aes(x = coast_dist/1000)) + # 海岸からの距離をkmで表示
  facet_wrap(~class, scales = "free_y", ncol = 1) +
  labs(x = "Distance from coast(km)")

標高差をヒストグラムで可視化

site %>%
  drop_na(paleo) %>% 
  ggplot() +
  geom_histogram(aes(x = Elv)) +
  facet_wrap(~class, scales = "free_y", ncol = 1) +
  scale_x_continuous(limits = c(0, 600)) +
  labs(x = "Elevation(m)")