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

# なぜRをつかうのか

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

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

# パッケージ読み込み

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

# データ読み込み

[北海道オープンデータポータル](https://www.harp.lg.jp/opendata/)の埋蔵文化財包蔵地データを利用します。ほかに行政界データを読み込みます。

```{r include=FALSE}
# データ読み込み
# 包蔵地ポイントデータ
site_pnt <- st_read("包蔵地_ポイント_20190913_JGD2011.shp")
# 包蔵地ポリゴンデータ
site_ply <- st_read("包蔵地_ポリゴン_20190913_JGD2011.shp")
# 市町村界
hkd <- st_read("N03-19_01_190101.shp")
```
この状態で地図を表示してみます。

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

![01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e160ef53-e408-a071-a3fd-08c46f26e8f9.png)

# ポリゴン重心を算出

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

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

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


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

# 厚沢部町域を抽出

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

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

![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/f36fb2ce-97b9-0af7-ac0c-26f03b8353fe.png)

# 座標系を変更

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

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

# 厚沢部町域の遺跡を抽出

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

```{r}
asbSite <- 
  st_intersection( x = asb3100 , y = site3100 )  # 町域（asb3100）で遺跡点(site3100)を切り抜き
# 地図を表示
asbSite %>%
  ggplot() +
    geom_sf(data = asb3100) +
    geom_sf() +
    theme_map() 
```
![03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/acd3386c-8091-19aa-807f-db4dbc9d7a69.png)

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

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

```{r}
site3100Na <-
  st_difference( x = site3100 , y = asb3100 ) 
# 地図を表示
site3100Na %>%
  ggplot() +
    geom_sf(data = hkd) +
    geom_sf() +
    theme_map()
```
厚沢部町域（赤いエリア）から遺跡が消えているのがわかります。
![05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/266f373f-2e7b-8023-cede-88c0088e449a.png)

# 時期ごとの遺跡分布図

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

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

```{r}
asbSite  %>%
  ggplot() +
    geom_sf(data = asb3100) +
    geom_sf( aes( colour = class ) , show.legend = "point") +
    theme_map() +
    scale_colour_viridis(option = "cividis" , discrete = TRUE)
```
![06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/5ffe8bb0-05fc-2cbc-1290-65801da4a528.png)

# 河川を読み込む

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

```{r}
wl <- 
  st_read("WL_polygon_utm54.shp" ,
    options = "ENCODING=CP932")   # 河川ポリゴン読み込み
wlasb <- 
  st_intersection( x = asb3100 , y = wl )  # 厚沢部町域で河川を切り抜き
# 地図を表示
wlasb %>%
  ggplot() +
    geom_sf()
```
![07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/08d25c0e-8b9b-3a16-fde6-a2a12076536d.png)

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

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

```{r}
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)
```
![08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/42e81fbf-aefe-d7ff-d36c-8936357bacfc.png)


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

```{r}
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"
    )
```
![09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/3d0ecde2-799b-588b-4136-da572e126f17.png)

# ボロノイポリゴンの作成

[こちら](https://www.jla-data.net/eng/creating-and-pruning-random-points-and-polygons/)を参考にしました。

```{r}
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で描画します。

```{r}
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() 
```

![10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e6589fe8-6d75-8b90-933e-88e27ca3fcb5.png)