考古学のためのRでGISーラスタ編

#  国土地理院基盤地図情報数値標高モデルをRで読み込む

`fgdr`パッケージを利用して、国土地理院基盤地図情報数値標高モデルを直接Rで読み込みます。

必要なパッケージ一式を読み込みます。

```{r}
library(tidyverse)
library(ggthemes)
library(sf)
library(viridis)
library(raster)
library(fgdr)
library(cptcity)
```
数値標高モデルをRのラスターオブジェクトに変換する手順は次のとおりです。

1. xmlファイルのリスト作成
2. `read_fgd_dem()`関数でxmlをrasterクラスに変換
3. for文でラスタオブジェクトを結合
4. Naが入力されているグリッド（海域＝-9999）を0に割り当て

# 基盤地図情報数値標高モデルを読み込み

以下のコードは、おおむね[Rで国土地理院 基盤地図情報データを扱う](https://uribo.hatenablog.com/entry/2018/12/16/084535)のコピペです。

```{r}
#1 ダウンロードして解答したxmlファイルのリストを作成します。
xml <- list.files(path="../01data/assabu_dem/", pattern=".*\\.xml", full.names=T) 
#2 
dem <- read_fgd_dem(xml[1], resolution = 10 , return_class = "raster")  # 最初のラスタ
#3 for文でxmlファイルをラスタオブジェクトに変換して結合
for(i in 2:length(xml) ) {
  dem1 <- read_fgd_dem(xml[i], 
               resolution = 10 ,
               return_class = "raster" )
  dem <- merge(dem , dem1,tolerance = 0.5)
}
#4  海域=-9999に0を代入
dem[dem<0] <- 0
# 地図を表示
plot(dem)
```
![01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/d2b6a5c1-b602-8651-8182-2f07103ed539.png)

# 座標系の変換

JGD2000UTMzone54に変換します。

```{r}
dem3100 <- 
  raster::projectRaster( dem , crs = CRS( "+init=epsg:3100" ) ) #JGD2000UTMzone54に変換
```

# 解像度を変更する

解像度を変更します。
10mメッシュのままでは、ggplotで描画する際に2分ほどかかってしまいます。
読み込んだDEMは9.21*12.30のメッシュなので、およそ50mメッシュに変更したいと思います。

```{r}
res(dem3100)  # 解像度を確認
[1]  9.21 12.30
```
`raster::aggregate()`関数でDEMの解像度を50mメッシュに変更します。`fact=5`でセルサイズを5倍にします。逆に`fact=-5`でセルサイズが5分1になります。

```{r}
dem3100_fif <-  
  dem3100 %>%
    raster::aggregate( fact=5 )  # fact=5 でセルを5倍にする
> res(dem3100_fif)  # 
[1] 46.05 61.50

```

# 切り抜き

厚沢部町域でラスタを切り抜きます。以下は前処理。

```{r}
hkd <- st_read("../01data/hkd_town_area/N03-19_01_190101.shp")   # 行政界データ読み込み
asb <- hkd %>% filter(N03_004 == "厚沢部町")    # 厚沢部町域を抽出
asb3100 <- asb %>%
  st_transform(3100)    # JGD2000 utm zon54  # 座標系をJGD2000UTMzone54に変換
```
厚沢部町域のbbox情報を抽出して、DEMのクロップ範囲に設定します。`raster::crop()`関数の引数で`raster::extent()`を指定し、クロップ範囲に厚沢部町域のbbox情報を指定します。町域の周囲に500mずつクリアランスを設定するため、minは500mを減じ、maxは500mを加えています。

```{r}
a <- st_bbox(asb3100)
     xmin      ymin      xmax      ymax 
 429609.0 4626320.7  456333.5 4656001.3 
dem3100_fif_cr <-
  dem3100_fif %>%
    raster::crop(   
      raster::extent( 
        c( a$xmin-500 , a$xmax+500 , a$ymin-500 , a$ymax+500  ) # 厚沢部町域のbbox情報をextent()に代入
      ) 
    ) 
  # extent()で指定する座標は順序に注意　
  # 500mずつクリアランスとる　
plot(dem3100_fif_cr)    # 地図を描画
```
![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/16e6e7cc-a9c7-44de-b4b9-49c3a3752e78.png)
   
# tibble形式に変換

ggplotで描画したいので、tibble形式に変換します。rasterオブジェクトからtibble形式には直接変換できないので、`as.data.frame()`関数でデータフレームに変換してからtibble形式に変換します（tbbleにしなくても描画できます）。

```{r}
df_dem <-
   dem3100_fif_cr %>%   # ラスターDEM
    as.data.frame( xy = TRUE ) %>% # rasterをデータフレームに変換
    tibble::as_tibble() %>% # tipple形式に変換
    dplyr::rename("Elevation" = layer)  # 標高値をElevationにリネーム
```

# ggplotで描画

標高ラスタ、等高線、河川、厚沢部町域、遺跡ポイントデータをggplotで描画します。カラーパレットには`cptcity`パッケージの`tp_tpushum`を使用しました。

```{r}
ggplot() +
  geom_raster(data = df_dem ,   # 標高ラスタ
            aes(x, y, fill = Elevation),
            hjust = 0, vjust = 0) +
  scale_fill_gradientn(colours = cptcity::cpt(pal = "tp_tpushum" , n = 50)) + # カラーパレットにcptを利用
  geom_contour(data = df_dem, aes(x, y, z = Elevation), col = "gray70", size = 0.2) +  # 等高線描画
  geom_sf(data = wlasb , fill = "#B4CDD8" , colour = "#B4CDD8") +  # 河川
  geom_sf(data = asb3100 , col = "gray40" , fill = NA, lwd = 0.8) +   # 厚沢部町域
  geom_sf(data=asbSite , aes( colour = class ) , show.legend = "point") +   # 遺跡
  theme_void() +
  scale_colour_viridis(option = "cividis" , discrete =TRUE ) +
  guides(fill=FALSE) 
```
分布図の邪魔にならないカラーパレットだと思います。
![03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e6e4ea60-487f-eed7-f158-cf657afbf933.png)

# rasterパッケージを駆使して地形指標を作成する

## 傾斜区分図
```{r}
slp <-
  dem3100_fif_cr %>%
    terrain( opt = "slope") 
plot(slp)   # 傾斜区分図を描画
```
![04slp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e4a83fbb-58cb-bae1-0c1e-39bb5a773597.png)

## 傾斜方位図
```{r}
asp <-
  dem3100_fif_cr %>%
    terrain( opt = "aspect") 
plot(asp)   # 傾斜方位図を描画
```
![05asp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/ddf8abe9-2676-7218-bb86-00390b56c05f.png)

## 陰影起伏図

```{r}
hill <-
  hillShade( slp , asp )
plot(hill)   # 陰影起伏図を描画
```
![06hill.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e221ca06-4d19-97de-4b5d-6866cfa1a2ed.png)

## TPI

Profilecurvatureとか地表開度と似た地形指標（だと思う）。

```{r}
# TPIは谷部で小さく尾根部で高い値をとる
tpi <-
  dem3100_fif_cr %>%
    terrain( opt = "TPI") 
plot(tpi)   # tpiを描画
```

![07tpi.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e3403ec5-dabb-b248-4b06-6898020caa98.png)

# 陰影起伏図を使った遺跡地図

## tibble形式に変換

```{r}
df_hill <-
   hill %>%   # ラスター陰影起伏図
    as.data.frame( xy = TRUE ) %>% # rasterをデータフレームに変換
    tibble::as_tibble() %>% # tipple形式に変換
    dplyr::rename("Shade" = layer)  # 標高値をElevationにリネーム
```

# ggplotで描画

標高ラスタに代えて、陰影起伏図を背景に使用します。

```{r}
ggplot() +
  geom_raster(data = df_hill ,   # 陰影起伏ラスタ
            aes(x, y, fill = Shade),
            hjust = 0, vjust = 0 , alpha = 0.5 ) +
  scale_fill_gradientn(colours=c("black","gray50","white"),limits=c(0,1)) +
  geom_contour(data = df_dem, aes(x, y, z = Elevation), col = "gray40", size = 0.2) +  # 等高線描画
  geom_sf(data = wlasb , fill = "#B4CDD8" , colour = "#B4CDD8") +  # 河川
  geom_sf(data = asb3100 , col = "gray40" , fill = NA, lwd = 0.8) +   # 厚沢部町域
  geom_sf(data=asbSite , aes( colour = class ) , show.legend = "point") +   # 遺跡
  theme_void() +
  scale_colour_viridis(option = "cividis" , discrete =TRUE ) +
  guides(fill=FALSE) 
```
![09hill.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e1f83bb6-29fc-bf7f-126d-440addf2fd9c.png)

# 傾斜区分図を使った遺跡地図

## tibble形式に変換

```{r}
df_slp <-
   slp %>%   # ラスター傾斜区分図
    as.data.frame( xy = TRUE ) %>% # rasterをデータフレームに変換
    tibble::as_tibble()  # tipple形式に変換
```

# ggplotで描画

```{r}
ggplot() +
  geom_raster(data = df_slp ,   # 標高ラスタ
            aes(x, y, fill = slope),
   ) +
  scale_fill_gradientn(colours = cptcity::cpt(pal ="colo_alpen_ski_slope_lunch"  , n = 50)) + # カラーパレットにcptを利用
  geom_sf(data = wlasb , fill = "#B4CDD8" , colour = "#B4CDD8") +  # 河川
  geom_sf(data = asb3100 , col = "gray40" , fill = NA, lwd = 0.8) +   # 厚沢部町域
  geom_sf(data=asbSite , aes( colour = class ) , show.legend = "point") +   # 遺跡
  theme_void() +
  scale_colour_viridis(option = "cividis" , discrete =TRUE ) +
  guides(fill=FALSE) 
```
![10slope.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/37f68e36-ff9b-d595-ec21-c246d87f1993.png)