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

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

fgdrパッケージを利用して、国土地理院基盤地図情報数値標高モデルを直接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で国土地理院 基盤地図情報データを扱うのコピペです。

#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)

座標系の変換

JGD2000UTMzone54に変換します。

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

解像度を変更する

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

res(dem3100)  # 解像度を確認
[1]  9.21 12.30

raster::aggregate()関数でDEMの解像度を50mメッシュに変更します。fact=5でセルサイズを5倍にします。逆にfact=-5でセルサイズが5分1になります。

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

切り抜き

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

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を加えています。

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)    # 地図を描画

tibble形式に変換

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

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を使用しました。

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) 

分布図の邪魔にならないカラーパレットだと思います。

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

傾斜区分図

slp <-
  dem3100_fif_cr %>%
    terrain( opt = "slope") 
plot(slp)   # 傾斜区分図を描画

傾斜方位図

asp <-
  dem3100_fif_cr %>%
    terrain( opt = "aspect") 
plot(asp)   # 傾斜方位図を描画

陰影起伏図

hill <-
  hillShade( slp , asp )
plot(hill)   # 陰影起伏図を描画

TPI

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

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

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

tibble形式に変換

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

ggplotで描画

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

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) 

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

tibble形式に変換

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

ggplotで描画

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)