Zenn
👏

R：sf 複数行のPOINTを1行のLINESTRINGに変換する、その逆も

2021/11/10に公開
R
GIS
gps
tech

はじめに

タイトルに書いたことをやろうとして、意外と手間取ったというかスマートにやる方法が見つからなかった。具体的にはGPSロガーで取得した点のデータを1つの線にする処理である。この記事には愚直な方法を備忘録として残す。
点（POINT）->線（LINESTRING）にするコードだけ先に書いておく。

# gps_sf_points は複数行にわたるsfオブジェクトで geometryはsfc_POINT
gps_sf_points %>% 
  sf::st_geometry() %>% 
  purrr::as_vector() %>% 
  matrix(ncol = 2, byrow = TRUE) %>% 
  sf::st_linestring() %>% 
  sf::st_sfc(crs = sf::st_crs(gps_sf_points)) # sfじゃなくてsfcだけどまあいいでしょう

GPSロガーのデータを読み込む

GPSロガーとしてAndroidのアプリ を使用した。このアプリは計測データをCSVとして出力できる。計測データは当然緯度経度を含むので、sfパッケージを利用して好きにいじれる。読み込みは次のように行う。出力をみると、そこそこな数の列を持つのがわかる。このうち、longitudeが経度で、latitudeが緯度である。

# 拡張子はtxtだけど中身はカンマ区切り
input <- "yyyyMMdd-hhmmss.txt"
gps_df <- readr::read_csv(input, col_names = T)
output
## ─ Column specification ──────────────────────────────────
## cols(
##   type = col_logical(),
##   `date time` = col_datetime(format = ""),
##   latitude = col_double(),
##   longitude = col_double(),
##   `accuracy(m)` = col_double(),
##   `altitude(m)` = col_double(),
##   `geoid_height(m)` = col_logical(),
##   `speed(m/s)` = col_double(),
##   `bearing(deg)` = col_double(),
##   sat_used = col_logical(),
##   sat_inview = col_double(),
##   name = col_character(),
##   desc = col_character()
## )

データフレームからsfオブジェクトに変換する。緯度経度を列として持つのでsf::st_as_sf()で素直に作れる。

# 経度、緯度の順で指定。''でくくる必要あり
gps_sf_points <- gps_df %>% 
  sf::st_as_sf(coords = c('longitude', 'latitude'), crs = 4326) # EPSG4326だよね？

データ数は次のように2943である。点->線->点で戻したとき数が一致するか確認のため出しておく。

nrow(gps_sf_points)
output
## [1] 2943

mapview::mapview()で地図に重ねてみる。移動の軌跡が点として描かれていることがわかる。

GPSロガーのデータを地図上に重ねたもの。点として認識されている。

点（POINT）から線（LINESTRING）へ

準備が終わったので本題の処理に入る。「はじめに」で示した処理を再掲する。ただし今度は変数に代入している。「関数ひとつでできるんじゃないの？」と思っていたのだがうまくいかず、結局ベクトルに変換してマトリックスを経由して…のような回りくどい処理になってしまっている。なお、コメントにも書いてあるようにgps_sfc_linestringはsfオブジェクトではなくてsfcオブジェクトである。

gps_sfc_linestring <- gps_sf_points %>% 
  sf::st_geometry() %>% 
  purrr::as_vector() %>% 
  matrix(ncol = 2, byrow = TRUE) %>% #byrowが必要
  sf::st_linestring() %>% 
  sf::st_sfc(crs = sf::st_crs(gps_sf_points)) # sfじゃなくてsfcだけどまあいいでしょう

1feature で、LINESTRINGになっているのがわかる

gps_sfc_linestring
output
## Geometry set for 1 feature 
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 137.2122 ymin: 36.70658 xmax: 137.2321 ymax: 36.72491
## Geodetic CRS:  WGS 84
## LINESTRING (137.2319 36.71248, 137.2319 36.7124...

地図に重ねれば確かに線としてつながっていることがわかる。

地図上に重ねたもの。きちんと線になっている

線（LINESTRING）から点（POINT）へ

今度は逆の処理をする。線にできたんだから逆に点に戻せるのでは？と思って試した。これもベクトルに変換してマトリックスを経由して…という処理になる。出力を見れば正しくPOINTで、2943行分のデータに戻っていることがわかる。これの最終結果はsfオブジェクトである。

gps_sfc_linestring %>% 
  sf::st_cast("MULTIPOINT") %>% 
  purrr::as_vector() %>% 
  matrix(ncol = 2) %>% # byrowは不要 
  data.frame() %>% 
  sf::st_as_sf(coords = c("X1", "X2"), crs = sf::st_crs(gps_sfc_linestring))
output
## Simple feature collection with 2943 features and 0 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 137.2122 ymin: 36.70658 xmax: 137.2321 ymax: 36.72491
## Geodetic CRS:  WGS 84
## First 10 features:
##                     geometry
## 1  POINT (137.2319 36.71248)
## 2  POINT (137.2319 36.71248)
## 3  POINT (137.2319 36.71248)
## 4  POINT (137.2319 36.71248)
## 5  POINT (137.2319 36.71248)
## 6  POINT (137.2319 36.71248)
## 7  POINT (137.2319 36.71248)
## 8  POINT (137.2319 36.71249)
## 9  POINT (137.2319 36.71249)
## 10  POINT (137.2319 36.7125)

まとめ

愚直な方法ではあるが、複数のPOINTとLINESTRINGとを相互変換する処理について整理した。

おまけ（失敗例）

素直にやるならこんな感じかな？と思って書いたコードはねらった結果にならなかった。

gps_sf_points %>% 
  dplyr::group_by(type) %>% # type列は全部TRUEなのでこれでグループ化 
  dplyr::summarise() %>% 
  sf::st_geometry() %>% 
  sf::st_cast("LINESTRING") %>% 
  mapview::mapview()


LINESTRINGではあるけど、もとの順番が崩れて経度順につながってしまっている

Discussion