Rのsfパッケージを利用して、ポリゴンの長軸と短軸を自動計測する

Rのsfパッケージを利用して、ポリゴンの長軸と短軸を自動計測します。

状況説明

次のような礎石配置図があり、この礎石の長軸と短軸を自動で計測します。

作戦

長軸と短軸を次のように定義します。

1. 長軸はポリゴンを構成する点群の分散が最大化するベクトル
2. 短軸は長軸と直行するベクトル

このように定義すると、長軸は主成分分析の主成分1のベクトル、短軸はこれと直交する主成分2のベクトルとして取り出せそうです。
そのうえで、各主成分スコアの最大値と最小値の差を取り出せば、長軸と短軸が求められそうです。

データ準備

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

# データ読込
data <-
  st_read("soseki.gpkg") %>% 
  filter(Layer == "10")    # Layerの10が外形線の場合
 
# データ表示
plot(data["Layer"]) 

ラインをポリゴンに変換する

今回のデータはラインデータで記録されていたのでこれをポリゴンに変換します。
st_polygonizeでポリゴン化するためにはラインの始点と終点が一致している必要があります。

# ラインをポリゴンに変換
data_poly <- st_polygonize(data)

ジオメトリを取り出す

対象となる礎石のジオメトリを取り出します。
st_geometryはジオメトリを抽出する関数です。

# 対象ジオメトリを抽出
geom <- st_geometry(data_poly[1, ])   # 対象礎石の1番目を取り出します

座標の抽出

st_coordinates関数でgeomの座標をmatrix形式で取り出します。

# 座標抽出
coords <- st_coordinates(geom)[, 1:2]

抽出したポリゴン各点の座標に対してPCA（主成分分析）を実行します。

# PCA実行
pca <- prcomp(
	coords, 
	center = TRUE,    # X Y座標を0を中心にシフト
	scale. = FALSE    # データの基準化は行わない
	)

長軸と短軸の長さを算出

# 各軸の範囲を算出 
axis_lengths <- 
	apply(pca$x,    # 各主成分スコア
		2,    # # 列ごとの処理 1だと行ごとの処理 
		function(x) diff(range(x))    # 最大値と最小値（range(x)）の差（diff()）を計算
		) 
# 長軸・短軸を表示
axis_lengths
# PC1が長軸、PC2が短軸
     PC1       PC2 
0.3687980 0.2830245 

長軸は0.368m、短軸は0.283mとなりました。

長軸短軸のベクトルとポリゴンオブジェクトを描画

実際に作成した長軸短軸の線分と、もとのポリゴンオブジェクトを重ねて描画します。

## 準備作業
# 中心点（重心）算出
center <- colMeans(coords)
# 主軸・副軸の単位ベクトルの抽出
major_axis_dir <- pca$rotation[, 1]  # 第1主成分（長軸方向）
minor_axis_dir <- pca$rotation[, 2]  # 第2主成分（短軸方向）

## ステップ1：長軸・短軸の端点を計算
# 長軸・短軸の中心点を基準に両端点を作成
long_axis_pts <- rbind(center + major_axis_dir * axis_lengths[1] / 2,
                       center - major_axis_dir * axis_lengths[1] / 2)

short_axis_pts <- rbind(center + minor_axis_dir * axis_lengths[2] / 2,
                        center - minor_axis_dir * axis_lengths[2] / 2)

## ステップ2：LINESTRINGに変換しsf化
# 長軸LINESTRING
long_axis_line <- st_linestring(long_axis_pts)
long_axis_sf <- st_sf(axis = "long", geometry = st_sfc(long_axis_line, crs = st_crs(data_poly_clean)))
# 短軸LINESTRING
short_axis_line <- st_linestring(short_axis_pts)
short_axis_sf <- st_sf(axis = "short", geometry = st_sfc(short_axis_line, crs = st_crs(data_poly_clean)))
# 両方結合
axis_sf <- rbind(long_axis_sf, short_axis_sf)

## ステップ3：ggplot2 で描画
ggplot() +
  geom_sf(data = data_poly_clean[3, ],
          fill = NA,
          color = "black") +  # ポリゴン1つだけ描画
  geom_sf(data = axis_sf, size = 1.2) +
  theme_classic()

重心点が少しずれているようですが、どこを計測したのかは読み取ることができます。
まあまあ、納得のいく結果です。

QGISで計測した結果

ぴったり一致するというわけではありませんが、自動計測した数値と比較してプラスマイナス1cm以内におさまっていますので、十分実用的と言えます。

すべてのデータをまとめて一括処理

すべての礎石データを一括して計測し、CSVファイルに書き出すコードです。

# データ読込
data <-
  st_read("soseki.gpkg") %>% 
  filter(Layer == "10")    # Layerの10が外形線の場合
# ラインをポリゴンに変換
data_poly <- st_polygonize(data)

# POLYGONまたはMULTIPOLYGONを含むものかどうかを判定する関数を定義
has_polygon <- function(geom) {
  # 空は除外
  if (st_is_empty(geom)) return(FALSE)
  
  # ジオメトリ型を取得
  gtype <- st_geometry_type(geom)
  
  # POLYGON or MULTIPOLYGON なら OK
  if (gtype %in% c("POLYGON", "MULTIPOLYGON")) return(TRUE)
  
  # GEOMETRYCOLLECTION なら POLYGON を抽出してチェック
  if (gtype == "GEOMETRYCOLLECTION") {
    poly_part <- st_collection_extract(geom, "POLYGON")
    return(length(poly_part) > 0)
  }
  
  # その他（POINT, LINESTRING など）は FALSE
  return(FALSE)
}

# ステップ2：filter で POLYGON を含む行だけ抽出
# TRUE のインデックスだけ抽出
poly_indices <- map_lgl(st_geometry(data_poly), has_polygon)

# フィルタリングして、ポリゴンのみを含む sf に
data_poly_clean <- data_poly[poly_indices, ]

# ステップ3：各ポリゴンの主成分分析を実行
## 空のデータフレーム作成
axis_lengths <- 
  data.frame(
    long = rep(NA_real_,
               nrow(data_poly_clean)), 
    short = rep(NA_real_,
                nrow(data_poly_clean))
  )

## 一括処理の実行
for(i in 1:nrow(data_poly_clean)){
  geom <- st_geometry(data_poly_clean[i, ])
  # GEOMETRYCOLLECTION の場合は最初の POLYGON を取り出す
  if (st_geometry_type(geom)[1] == "GEOMETRYCOLLECTION") {
    geom_parts <- st_collection_extract(geom, "POLYGON")
    if (length(geom_parts) == 0) {
      stop("GEOMETRYCOLLECTION に POLYGON が含まれていません")
    }
    geom <- geom_parts[[1]]
  }
  # 座標抽出
  coords <- st_coordinates(geom)[, 1:2]
  # PCA実行
  pca <- prcomp(coords, center = TRUE, scale. = FALSE)
  # 中心点（重心）
  center <- colMeans(coords)
  
  # 主軸・副軸の単位ベクトル
  major_axis_dir <- pca$rotation[, 1]  # 第1主成分（長軸方向）
  minor_axis_dir <- pca$rotation[, 2]  # 第2主成分（短軸方向）
  
  # 各軸の標準偏差（scoreの散らばり） → 全長にするには×2
  axis_lengths[i,] <- 
    apply(pca$x,    # 各主成分スコア
          MARGIN = 2,    # 列ごとの処理 1だと行ごとの処理
          function(x) diff(range(x)))    # 最大値と最小値（range(x)）の差（diff()）を計算
}

# axsis_lengths の列名を設定
axis_lengths_No <-
  axis_lengths %>% 
  as_tibble() %>%
  bind_cols(data_poly_clean$No) %>% 
  rename(No = ...3) %>%
  arrange(No) %>% 
  mutate(
    long = long * 100,  # 長軸の長さをcmに変換
    short = short * 100 # 短軸の長さをcmに変換
  ) %>%
  mutate(
    long = round(long, 1), 小数点一桁まで表示
    short = round(short, 1)
  )

# csvに出力
write_csv(axis_lengths_No, 
          "axis_lengths.csv")

このようなデータフレームが生成され、長軸と短軸を取り出すことに成功しました。
251個の礎石のうち、始点と終点が一致する（ポリゴン化できる）230個の礎石の長軸と短軸の長さを自動的に計測できています。

axis_lengths_No
# A tibble: 230 × 3
    long short No   
   <dbl> <dbl> <chr>
 1  33    31.8 001  
 2  29.5  26.4 002  
 3  36.9  28.3 003  
 4  33.3  23.1 006  
 5  42.2  22   007  
 6  46.5  29.6 008  
 7  44.9  45.4 009  
 8  42.9  31.8 010  
 9  34.8  34.1 011  
10  35.9  22.3 012  
# ℹ 220 more rows
# ℹ Use `print(n = ...)` to see more rows