Rとterraパッケージを使った断面図の作成

地理情報システム（GIS）を利用して地形の断面図を作成することは、地形解析や土木設計において非常に重要です。QGISのProfile toolは断面図を作成する優れたプラグインですが、出力されるのは画像形式のみです。これをCADソフトウェアなどで扱えるようにするために、Rのterraパッケージを用いてCSV形式の断面を作成する方法を紹介します。

# パッケージ読み込み
library(terra)
library(tidyterra)
library(tidyverse)
library(ggthemes)
library(cptcity)

データの準備

使用する地形データはGeoTIFF形式のデジタル標高モデル（DEM）です。また、断面を切るためのラインデータはGeoPackage形式で用意します。

使用するのはこのような地形データで青いラインに沿って断面を切ります。

# GeoTIFFファイルの読み込み
dem <- rast("DEM.tif")
# セクションラインデータの読み込み
line <- vect("SectionLine.gpkg")

サンプリングと標高の抽出

ライン上を等間隔にサンプリングし、各点の標高を抽出します。

# ライン上を等間隔にサンプリング
sampled_points <-
  line %>% 
  spatSample(size = 100, method = "regular")    # サンプリングサイズ100点
# 標高を抽出
profile_values <- extract(dem, sampled_points)

距離計算

サンプリングした点の座標を取得し、原点からの距離を計算します。

# 距離計算
## サンプリングした点の座標を取得
coords <- crds(sampled_points) 
# 距離を計算するために、サンプリングした点の座標を取得
distances <- c(0, cumsum(sqrt(rowSums(diff(coords)^2))))

データの保存

計算したデータをデータフレームに格納し、CSVファイルとして保存します。

# 距離計算したデータをtibbleに格納
profile_df <- tibble(
  x = coords[,1],
  y = coords[,2],
  distance = distances,
  elevation = profile_values[[2]]
)

# データフレームをCSVファイルに保存
write_csv(profile_df, "profile.csv")

データの確認

profile_dfはx座標、y座標、原点からの距離、各地点の標高値で構成されています。

> profile_df
# A tibble: 100 × 4
       x        y distance elevation
   <dbl>    <dbl>    <dbl>     <dbl>
 1 7657. -235788.     0         78.4
 2 7658. -235783.     4.89      78.1
 3 7658. -235778.     9.77      78.1
 4 7658. -235773.    14.7       78.1
 5 7658. -235768.    19.5       78.0
 6 7658. -235764.    24.4       77.9
 7 7658. -235759.    29.3       77.7
 8 7658. -235754.    34.2       77.6
 9 7659. -235749.    39.1       77.3
10 7659. -235744.    44.0       77.1
# ℹ 90 more rows
# ℹ Use `print(n = ...)` to see more rows

プロットの作成

最後に、ggplot2を用いて断面図をプロットし、視覚的に確認します。

# ggplotで断面をプロット
profile_df %>%
  ggplot(aes(x = distance, y = elevation)) +
  geom_line() +
  labs(title = "Elevation Profile",
       x = "Distance (m)",
       y = "Elevation (m)") +
  theme_bw() +
  coord_fixed(ratio = 2)    # Y軸を2倍に強調する