考古学のためのRでGIS-ベクタ編

なぜRをつかうのか

今回紹介する方法は、QGISならとても簡単です。QGISのようなGISソフトを使えば簡単にできることをあえてRで行う理由は、研究の再現性です。

単純にドットを表示するだけならともかく、データの前処理から分析を経て出図にいたるプロセスは、その手順が長くなればなるほど、再現できなくなっていきます。Rを使うことで、複雑なプロセスを根気よくたどり、研究を再現することが可能になります。

パッケージ読み込み

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

データ読み込み

北海道オープンデータポータルの埋蔵文化財包蔵地データを利用します。ほかに行政界データを読み込みます。

# データ読み込み
# 包蔵地ポイントデータ
site_pnt <- st_read("包蔵地_ポイント_20190913_JGD2011.shp")
# 包蔵地ポリゴンデータ
site_ply <- st_read("包蔵地_ポリゴン_20190913_JGD2011.shp")
# 市町村界
hkd <- st_read("N03-19_01_190101.shp")

この状態で地図を表示してみます。

# 地図を表示
ggplot() +
  geom_sf(data = hkd ) +
  geom_sf(data = site_pnt) +
  geom_sf(data = site_ply ) +
  theme_map()

ポリゴン重心を算出

遺跡ポリゴンデータとポイントデータを一つのオブジェクトとして扱いために、ポリゴンデータをポイントデータに変換します。st_centroid()関数で重心点を算出します。

# st_centroid()関数でポリゴン重心点を算出
site_plyToPnt <- site_ply %>% 
  sf::st_centroid() # ポリゴン重心を算出

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

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

厚沢部町域を抽出

行政界データから厚沢部町域のみを抽出します。dplyr::filter()で厚沢部町を抽出します。

asb <- hkd %>% filter(N03_004 == "厚沢部町")
# 地図を表示
ggplot() +
  geom_sf(data = asb ) +
  theme_map()

座標系を変更

st_transform()関数で緯度経度系から投影座標系に変換します。
緯度経度系のままでは、空間結合系の機能が使えません。
EPSGコードを引数に指定できるので、JGD2000 utm zon54に該当する3100を指定します。

site3100 <- site %>% 
  st_transform(3100)    # JGD2000 utm zon54
asb3100 <- asb %>%
  st_transform(3100)    # JGD2000 utm zon54  

厚沢部町域の遺跡を抽出

st_intersection()は「交差」に該当する空間操作です。
第1引数に切り取る範囲を示すポリゴン（厚沢部町域=asb3100）を指定し、第2引数に切り取られる側の地物（遺跡点=site3100）を指定します。

asbSite <- 
  st_intersection( x = asb3100 , y = site3100 )  # 町域（asb3100）で遺跡点(site3100)を切り抜き
# 地図を表示
asbSite %>%
  ggplot() +
    geom_sf(data = asb3100) +
    geom_sf() +
    theme_map() 

厚沢部町域以外の遺跡を抽出

st_difference()はst_intersection()の補集合にあたります。
厚沢部町域以外の遺跡点が抽出されます。
第1引数に切り取られる側の地物（遺跡点）を指定し、第2引数に切り取る範囲を示すポリゴン（厚沢部町域）を指定します。
ポリゴンとポイントの引数がst_intersection()と逆になるので注意が必要です。

site3100Na <-
  st_difference( x = site3100 , y = asb3100 ) 
# 地図を表示
site3100Na %>%
  ggplot() +
    geom_sf(data = hkd) +
    geom_sf() +
    theme_map()

厚沢部町域（赤いエリア）から遺跡が消えているのがわかります。

時期ごとの遺跡分布図

遺跡点のasbSite$classは遺跡の時期を示すデータ（筆者作成。道庁の包蔵地データには含まれていない）です。早期からアイヌ文化期までの8期に分類しています。

levels(asbSite$class)
[1] "早期"         "前期"         "中期"         "後期"         "晩期"        
[6] "続縄文"       "擦文文化期"   "アイヌ文化期"

asbSite  %>%
  ggplot() +
    geom_sf(data = asb3100) +
    geom_sf( aes( colour = class ) , show.legend = "point") +
    theme_map() +
    scale_colour_viridis(option = "cividis" , discrete = TRUE)

河川を読み込む

基盤地図情報から作成した河川ポリゴンを読み込み、st_intersection()を使って厚沢部町域ポリゴンで切り抜きます。

wl <- 
  st_read("WL_polygon_utm54.shp" ,
    options = "ENCODING=CP932")   # 河川ポリゴン読み込み
wlasb <- 
  st_intersection( x = asb3100 , y = wl )  # 厚沢部町域で河川を切り抜き
# 地図を表示
wlasb %>%
  ggplot() +
    geom_sf()

遺跡、町域、河川を重ねる

遺跡、町域、河川のデータを重ねて表示します。
遺跡は時期ごとに色を分けてみます。

asbSite %>%
  ggplot() +
    geom_sf(data = asb3100 , fill = NA) +　# 町域ポリゴン　塗りつぶしなし
    geom_sf(data = wlasb  , fill = "Grey" , colour = "Grey") +　# 河川ポリゴン　
    geom_sf( aes( colour = class ) , show.legend = "point") +　# 遺跡点データ　show.legend = "point"で凡例を点にする
    theme_map() +
    scale_colour_viridis(option = "viridis" , discrete = TRUE)

時期ごとにファセットします。

asbSite  %>%
  ggplot() +
    geom_sf(data = asb3100 , fill = NA) +
    geom_sf(data = wlasb  , fill = "Grey") +
    geom_sf( aes( colour = class ) , show.legend = "point") +
    facet_wrap(~class) +
    theme_map() +
    scale_colour_viridis(option = "viridis" , discrete = TRUE) +
    theme(
      legend.position= "none"
    )

ボロノイポリゴンの作成

こちらを参考にしました。

vor  <- asbSite %>%  # 遺跡地点データ
  st_geometry() %>% # 属性を全部捨ててジオメトリだけにする
  st_union() %>% # 行を全部結合
  st_voronoi() %>% # ボロノイを作成
  st_collection_extract(type = "POLYGON") %>% # オブジェクトのタイプを選択 type = "POLYGON"
  st_sf(crs = 3100) %>% # CRSを設定　JGD2000 utm zone54（crs=3100） を指定
  st_intersection(asb3100) %>% # 厚沢部町域でクリップ
  st_join(asbSite) %>% # もともとの遺跡地点のデータと結合（捨てた属性が元に戻る）
  st_set_agr("aggregate") # clean up

ggplotで描画します。

ggplot() +
  geom_sf(data = wlasb  , fill = "gray80" , colour = "gray80" ) +
  geom_sf(data = vor, col = 'gray60', fill = NA) +
  geom_sf(data = asb3100, col = 'gray25', fill = NA, lwd = 1) +
  geom_sf(data = asbSite, show.legend = "point" ,
    aes(colour = class) ) +
    scale_colour_viridis(option = "viridis" , discrete = TRUE) +
  theme_map()