计算R中两个sf点的夹角来计算道路线串的方向答案

【问题标题】：Calculate the angle between two sf points in R to calculate the direction of road line strings计算R中两个sf点的夹角来计算道路线串的方向
【发布时间】：2022-08-10 18:15:18
【问题描述】：

尝试使用每个 sf 线串的起点和终点来计算线串的 sf 向量的角度。我不想拆线。

我在这个类似的例子中使用了公式来计算度数 Calculate angle between two Latitude/Longitude points

下面的示例下载了一个小型道路网络，创建了一个函数 \'angles\'，它在数据框中添加一个带有角度的列，然后在 mapview 中绘制它，因此每条线都按角度着色。但是，这些值不正确。我认为公式的简单错误？

library(osmdata)
library(sf)
library(dplyr)


bb <- st_bbox(c(xmin = 5.22, xmax = 5.246, ymax = 52.237, ymin = 52.227), crs = st_crs(4326))

x <- opq(bbox = bb) %>%
  add_osm_feature(key = c(\'highway\')) %>%
  osmdata_sf()

##extract building polygons
roads <- x$osm_lines %>% 
  filter(highway == \"motorway\") %>% 
  select(osm_id, geometry)

angles <- function(x){
  
  ## define line finish coordinates
  f <- x %>%
    st_transform(28992) %>% 
    st_line_sample(sample = 1) %>%
    st_cast(\"POINT\") %>% 
    st_transform(4326) %>% 
    st_coordinates()
  
  ## define line start coordinates
  s <- x %>%
    st_transform(28992) %>% 
    st_line_sample(sample = 0) %>% 
    st_cast(\"POINT\") %>% 
    st_transform(4326) %>% 
    st_coordinates()
  
  ## get latitude of finish points
  lat2 <- f[,2]
  ## get latitudes of start points
  lat1 <- s[,2]
  

  
  ## get longitudes of start points
  lon2 <- f[,1]
  ## get longitudes of start points
  lon1 <- s[,1]
  
  ## delta longitudes
  #dlon <- f[,1]-s[,1]
  #theta <- (atan2(sin(dlon)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon)))*180/pi
  
  theta <- atan2(lat2-lat1, lon2-lon1)*180/pi
  
  x$angle_deg <- theta
  
  return(x)
  
}

a <- angles(roads) %>% 
  select(angle_deg)

library(mapview)

mapview(a)

标签： r sf angle

【解决方案1】：

我建议您考虑lwgeom::st_geod_azimuth() - 它会提供方位（以弧度为单位）。

需要考虑两件事：

它需要几何类型POINT - 因为从一个点到另一个点的方位是非常清晰的，但是LINESTRING 的斜率可能/将会改变一个部分到下一个部分
它将给出一个比原始点少一个元素的向量 - 最后一个点没有定义的方位，因为没有“下一个”点；为了解决这个问题，我在最后一个元素中添加了人工NA。

如果您希望以度为单位输出，只需将计算传递给units::set_units("degrees") 调用。

library(osmdata)
library(sf)
library(dplyr)


bb <- st_bbox(c(xmin = 5.22, xmax = 5.246, ymax = 52.237, ymin = 52.227), crs = st_crs(4326))

x <- opq(bbox = bb) %>%
  add_osm_feature(key = c('highway')) %>%
  osmdata_sf()

##extract building polygons
roads <- x$osm_lines %>% 
  filter(highway == "motorway") %>% 
  select(osm_id, geometry)  %>% 
  st_cast("POINT") %>% 
  mutate(bearing = c(lwgeom::st_geod_azimuth(.),
                     units::set_units(NA, "radians"))) # 1 extra point align vector length

mapview::mapview(roads["bearing"])

【讨论】：

非常感谢。我意识到也许我对这个问题并不完全清楚。我不想拆分线串，因为我正在使用角度来使用来自其他几何图形的数据填充这些线。但是，您的代码帮助很大。这是我拥有的最终代码，它可以调整您的代码以产生每个唯一线串的角度
@B_K oki，很高兴为您服务！

【解决方案2】：

非常感谢。我意识到也许我对这个问题并不完全清楚。我不想拆分线串，因为我正在使用角度来使用来自其他几何图形的数据填充这些线。但是，您的代码帮助很大。这是我拥有的最终代码，它可以调整您的代码以产生每个唯一线串的角度。

library(osmdata)
library(sf)
library(dplyr)

    bb <- st_bbox(c(xmin = 5.22, xmax = 5.246, ymax = 52.237, ymin = 52.227), crs = st_crs(4326))
    
    x <- opq(bbox = bb) %>%
      add_osm_feature(key = c('highway')) %>%
      osmdata_sf()
    
pts <- x$osm_lines %>% ## extract polylines from osm data
  filter(highway == "motorway") %>% ## only motorway network
  select(osm_id, geometry)  %>% ## remove excess columns
  st_cast("POINT") %>% ## split into constituent points
  group_by(osm_id) %>% 
  filter(row_number()==1 | row_number()==n()) %>% ## take start and end points of each osm id
  mutate(bearing = c(lwgeom::st_geod_azimuth(geometry))) %>% ## calculate bearing of start and end points
  mutate(bearing = units::set_units(bearing, "degrees")) %>% ## convert to degrees
  distinct(osm_id, .keep_all = TRUE) ## just take one point for each osm id
    
    roads <- x$osm_lines %>% 
      filter(highway == "motorway") %>% 
      select(osm_id, geometry)  %>% 
      mutate(bearing = pts$bearing)
    
    mapview::mapview(roads["bearing"])

【讨论】：

【解决方案3】：

您可能有兴趣将sfnetworks 包用于空间网络。它包括函数edge_azimuth() 来计算网络中边缘的方位角。

请注意，OSM 数据在形成干净、可路由的网络之前可能需要进行一些预处理。参见例如sfnetworks 关于预处理和清洁的小插曲：https://luukvdmeer.github.io/sfnetworks/articles/sfn02_preprocess_clean.html

代表：

library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(sf)
#> Linking to GEOS 3.10.1, GDAL 3.4.0, PROJ 8.2.0; sf_use_s2() is TRUE
library(dplyr)
library(sfnetworks)

bb <- st_bbox(c(xmin = 5.22, xmax = 5.246, ymax = 52.237, ymin = 52.227), crs = st_crs(4326))

x <- opq(bbox = bb) %>%
  add_osm_feature(key = c('highway')) %>%
  osmdata_sf()

net <- x$osm_lines %>% 
  filter(highway == "motorway") %>% 
  as_sfnetwork() %>%
  activate("edges") %>%
  mutate(bearing = edge_azimuth()) %>%
  mutate(bearing = units::set_units(bearing, "degrees"))

plot(st_as_sf(net)[, "bearing"])

^{由reprex package (v2.0.1) 创建于 2022-08-10}

【讨论】：