rcarbonパッケージは放射性炭素年代測定のデータベースを分析するためのパッケージです。較正年代を提供するだけではなく、様々な統計的検定にも対応しています。この記事では、rcorbonパッケージの空間分析機能について解説します。
全コード
北海道のC14年代の空間分析を行うための全コードです。
library(sf)
library(tidyverse)
library(spatstat)
library(rcarbon)
# read data
site <- st_read("data/HKDsiteMerge_JGD2000utm54.gpkg")
radio <- read_csv("data/C14HKD.csv")
# data merge
site_radio <-
left_join(site, radio, by = c("名称" = "遺跡名")) %>%
dplyr::select(名称, 種別, C14年代, C14年代シグマ, geom) %>%
rename(site = 名称, type = 種別, C14Age = C14年代, C14SD = C14年代シグマ)
# coordinate
coor <- st_coordinates(site_radio$geom)
# bind
site_radio_df <-
bind_cols(site_radio, coor) %>%
as.data.frame() %>%
dplyr::select(-geom) %>%
drop_na(C14Age, C14SD)
# Load the window polygon
hkd <- st_read("data/HokkaidoPolygon_JGD2000utm54.gpkg")
# Convert to spatstat owin object
hkd_owin <- as.owin(hkd)
## Calibrate and bin
x <-
calibrate(x = site_radio_df$C14Age,
errors = site_radio_df$C14SD,
normalised = FALSE,
verbose = FALSE)
bins1 <- binPrep(sites = site_radio_df$site,
ages = site_radio_df$C14Age,
h = 50) # 指定した数値の数でデータを分割
## Create centennial timeslices (see help doc for further argument details)
stkde1 <- stkde(x = x,
coords = site_radio_df[,c("X", "Y")],
win = hkd_owin,
sbw = 40000, # ガウシアンカーネルの標準偏差 たぶん空間距離を示す？
cellres = 5000, # グリッドの解像度 1ピクセルの距離
focalyears = seq(50000, 50, -100), # 時間スライスの指定
tbw = 50, # ガウシアンカーネルの標準偏差 たぶん時間軸の距離を示す？
bins = bins1,
backsight = 200, # 過去との比較タイムステップを指定する
outdir = tempdir(),
amount = 1,
verbose = FALSE)
# plot
par(mar = c(0.5, 0.5, 2.5, 2))
plot(x = stkde1,
focalyears = 500, # 対象とする時間スライスの指定 calBP
type = "all")
データについて
- 国立歴史民俗博物館が公開している遺跡調査報告書放射性炭素年代測定データベース
- 遺跡位置情報は北海道庁が公開している埋蔵文化財包蔵地（GISデータ）【北海道】
- 分析対象領域を定義する北海道のポリゴンデータは国土数値情報の行政区域北海道版
パッケージ読み込み
sf、tidyverse、spatstat、rcarbonの各パッケージを読み込みます。
spatstatは分析対象領域を示すowinデータを作成するために使用します。
# laod package
library(sf)
library(tidyverse)
library(spatstat)
library(rcarbon)
データの読み込み
# read data
site <-
st_read("data/HKDsiteMerge_JGD2000utm54.gpkg") # 北海道の遺跡GISデータ
radio <-
read_csv("data/C14HKD.csv") # C14年代測定データ
siteオブジェクトに格納したのは北海道の遺跡のGISデータで、radioオブジェクトに格納したのは遺跡調査報告書放射性炭素年代測定データベースからダウンロードした北海道のC14年代測定結果一覧表です。
siteの中身はこんな感じです。
site
Simple feature collection with 12423 features and 8 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 369758.1 ymin: 4584095 xmax: 889918.5 ymax: 5040627
Projected CRS: JGD2000 / UTM zone 54N
First 10 features:
名称 登載番号 所在地 種別 時代 立地 標高 出土 geom
1 坊主山遺跡 A-02-015 江別市対雁1・2 墳墓 続縄文,擦文,アイヌ 石狩川左岸段丘 10 m 土器、鉄器 POINT (543667.1 4774414)
2 江別第2チャシ跡 A-02-105 江別市対雁125-1ほか チャシ跡 アイヌ 旧豊平川右岸段丘端部 20 m POINT (542825.3 4774351)
3 美原遺跡 A-02-034 江別市美原河川敷地 遺物包含地 不明 POINT (546291.3 4774576)
4 元江別9遺跡 A-02-086 江別市元江別本町35･37 遺物包含地 縄文(中期) POINT (542836.3 4772974)
5 元江別5遺跡 A-02-059 江別市元江別810･874 遺物包含地 縄文 POINT (542700.6 4773580)
6 元江別4遺跡 A-02-058 江別市元江別 遺物包含地 縄文 POINT (542688.6 4773675)
7 元江別3遺跡 A-02-057 江別市元江別806 集落跡 縄文(晩期),続縄文 POINT (542722.3 4773775)
8 町村農場4遺跡 A-02-056 江別市元江別775･806 遺物包含地 縄文(中期) POINT (542866.7 4773846)
9 町村農場3遺跡 A-02-055 江別市対雁111～113 遺物包含地 縄文 POINT (543006.3 4774015)
10 対雁遺跡 A-02-052 江別市対雁97,100 遺物包含地 縄文(中期),続縄文 POINT (543296.5 4774423)
radioの中身はこんな感じ。3,792件の年代測定データが格納されています。
radio
# A tibble: 3,702 × 29
項番 番号 都道府県 遺跡名 所在地 サンプル採取地点等 試料の種類 試料コード1 試料コード2 時代 時代詳細 時代コード 試料番号 測定方法 C14年代 C14年代シグマ
<dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
1 1 33124 01：北海道 C143遺跡 札幌市中央区… 7a層 C143-02 土層（黒… T k 縄文… ---- B Beta-16… AMS法 3150 40
2 2 33125 01：北海道 C143遺跡 札幌市中央区… 9a層 C143-03 土層（黒… T k 縄文… ---- B Beta-16… AMS法 3900 40
3 3 33123 01：北海道 C143遺跡 札幌市中央区… 5層 出土遺物なし … 土層（黒… T k 続縄… ---- I Beta-16… AMS法 2240 40
4 4 33127 01：北海道 C424遺跡A地点 札幌市中央区… 2～3a層、竪穴住居… 炭化材 T f 擦文… ---- J Beta-13… AMS法 1290 40
5 5 33132 01：北海道 C424遺跡A地点 札幌市中央区… 7層 04-04区、埋没… クルミ T e 縄文… ---- B Beta-13… AMS法 2450 40
6 6 33131 01：北海道 C424遺跡A地点 札幌市中央区… 2～3a層、竪穴住居… 炭化材 T f 擦文… ---- J Beta-13… β線法 1260 70
7 7 33130 01：北海道 C424遺跡A地点 札幌市中央区… 2～3a層、竪穴住居… クルミ T g 擦文… ---- J Beta-13… AMS法 970 40
8 8 33129 01：北海道 C424遺跡A地点 札幌市中央区… 2～3a層 竪穴住居… 炭化材 T f 擦文… ---- J Beta-13… β線法 900 80
9 9 33128 01：北海道 C424遺跡A地点 札幌市中央区… 2～3a層 竪穴住居… 炭化材 T f 擦文… ---- J Beta-13… AMS法 880 40
10 10 33126 01：北海道 C424遺跡A地点 札幌市中央区… 2～3a層、竪穴住居… 炭化材 T f 擦文… ---- J Beta-13… AMS法 1140 40
# ℹ 3,692 more rows
# ℹ 13 more variables: 暦年較正用C14年代 <chr>, `暦年較正用C14年代±` <chr>, `δ13C（AMS）` <chr>, `δ13C（AMS）±` <chr>, `δ13C（IRMS）` <chr>, `分析者（著者）` <chr>,
# 分析機関 <chr>, 刊行年 <chr>, 報告タイトル <chr>, ページ <chr>, 備考 <chr>, 報告書名 <chr>, 発行者 <chr>
# ℹ Use `print(n = ...)` to see more rows
遺跡データと年代測定データを結合する
siteに格納した北海道の遺跡データとradioに格納した年代測定データを結合します。結合はleft_join関数を用い、遺跡名をキーにして強引に結合します。遺跡名が異なっている場合（「〇〇遺跡〇〇地点となっているケースなど）は結合できません。
# data merge
site_radio <-
left_join(site, radio, by = c("名称" = "遺跡名")) %>%
dplyr::select(名称, 種別, C14年代, C14年代シグマ, geom) %>%
rename(site = 名称, type = 種別, C14Age = C14年代, C14SD = C14年代シグマ)
次に、st_coordinates関数を用いて、遺跡位置情報をX、Y座標として抽出します。
# XY座標を抽出する
coor <- st_coordinates(site_radio$geom)
# XY座標だけが抽出される。
head(coor)
X Y
[1,] 542061 4775314
[2,] 542061 4775314
[3,] 542061 4775314
[4,] 542061 4775314
[5,] 542061 4775314
[6,] 542061 4775314
- 年代測定データと遺跡のXY座標を結合します。
- rcorbonではtibble形式でエラーが出ることがあったので、データフレーム形式にします。
- 最後にNA行を除去します。
# 年代測定データとXY座標を結合する
site_radio_df <-
bind_cols(site_radio, coor) %>% # 年代測定・遺跡データとXY座標を結合
as.data.frame() %>% # データフレーム形式に変換
dplyr::select(-geom) %>% # geometoryはいらないので消す
drop_na(C14Age, C14SD) # C14年代と標準偏差のNA行を除外する
# 最終的なデータ
head(site_radio_df)
site type C14Age C14SD X Y
1 対雁2遺跡 遺物包含地 2475 20 542061 4775314
2 対雁2遺跡 遺物包含地 2485 20 542061 4775314
3 対雁2遺跡 遺物包含地 2460 20 542061 4775314
4 対雁2遺跡 遺物包含地 2465 20 542061 4775314
5 対雁2遺跡 遺物包含地 2440 20 542061 4775314
6 対雁2遺跡 遺物包含地 2405 20 542061 4775314
owinオブジェクトを作成する
分析対象領域を定義するowinオブジェクトを作成します。国土数値情報の市町村界のGISデータを加工して北海道全域のポリゴンデータを作成しておきます。
# 北海道のポリゴンデータを読み込む
hkd <- st_read("data/HokkaidoPolygon_JGD2000utm54.gpkg")
# owin形式に変換する
hkd_owin <- as.owin(hkd)
plot(hkd_owin) # owinオブジェクトをプロット
作成したowinオブジェクトはこんな感じです。
年代較正
rcarbon::calibrate関数は年代較正を行います。測定年代と標準誤差を引数にとります。
x <-
rcarbon::calibrate(x = site_radio_df$C14Age, # C14年代
errors = site_radio_df$C14SD, # C14年代の標準偏差
normalised = FALSE,
verbose = FALSE)
ビニング
ビニング（binning）はちょっとわかりにくい概念ですが、一つの遺跡で飛び抜けて測定結果が多い場合に計測値をグループ化する機能です。突出した計測データ量をもつ遺跡の影響を回避するために用いられるようです。内部的にクラスター分析をしているようで、hオプションでグループ化の度合いを制御しています。
bins1 <-
rcarbon::binPrep(sites = site_radio_df$site,
ages = site_radio_df$C14Age,
h = 50) # クラスター分析のグループ化の度合いを制御する（らしい）
タイムスライスを作成する
いよいよ空間分析に向けたタイムスライスを作成します。タイムスライスとは、特定の年代の計測値の分布密度のことです。C14年代は時間的にも広がりをもつ概念であるため、単に空間補完を行えば良いわけではありません。時空間の2軸にわたって広がりと強度を評価するための計算を実行するのがrcarbon::stkde関数です。
この関数が優れているのは、単に特定の時点の空間分布を示すだけではなく、分布全体における特異性や1スパン前との比較も同時に行ってくれることです。
stkde1 <-
rcarbon::stkde(x = x,
coords = site_radio_df[,c("X", "Y")], XY座標を指定
win = hkd_owin,
sbw = 40000, # ガウシアンカーネルの標準偏差 たぶん空間距離
cellres = 5000, # グリッドの解像度 1ピクセルの距離
focalyears = seq(30000, 0, -250), # 時間スライスの指定 3万年前から現在までを250年ステップでスライス
tbw = 50, # ガウシアンカーネルの標準偏差 たぶん時間軸の距離
bins = bins1,
backsight = 500, # 過去との比較タイムステップを指定す
outdir = tempdir(),
amount = 1,
verbose = FALSE)
計算時間はけっこう長いです。
描画する
3000年前（縄文時代後期末）の遺跡分布状況です。
par(mar = c(0.5, 0.5, 2.5, 2))
plot(x = stkde1,
focalyears = 3000, # 対象とする時間スライスの指定 calBP
type = "all")
左から当該時期の分布密度、全体の分布密度、全体密度あたりの当該期の相対分布密度、500年前と比較した密度の増減です。時期ごとに密度図を作成することで、遺跡動態が見えてきそうです。
