R：sf 複数行のPOINTを1行のLINESTRINGに変換する、その逆も
はじめに
タイトルに書いたことをやろうとして、意外と手間取ったというかスマートにやる方法が見つからなかった。具体的には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)
## ─ 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)
## [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
## 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))
## 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