计算R中两个SF点之间的角度以计算路线线的方向

Question

试图使用每个SF Linestring的起点和终点来计算线条载体的SF矢量角度。我不想将阵容分开。

我在此类似示例中使用了该公式进行学位计算 计算两个纬度/经度点之间的

角度网络，创建一个函数“角度”，该函数将列添加到数据框架中，然后将其绘制在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)

Answer 1

我建议您考虑 lwgeom :: st_geod_azimuth（） - 它将给出轴承（以弧度为单位）。

要考虑的两件事：

• 它需要几何类型点 - 因为从一个点到另一个点是清晰的，但是 linestring的斜率 may/ther将更改一个部分
• 它将给出一个元素的向量小于原始点 - 最后一点没有定义的轴承，因为没有“下一个”点；为了解决这个问题，我将人造 na 添加到最后一个元素。

如果您希望以学位的输出，只需将计算输出到a units :: set_units :: set_units :: set_units（“ demrees”） 呼叫。

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"])

Answer 2

非常感谢。我意识到也许我对这个问题并不完全清楚。我不想在使用角度将其与其他几何形状的数据填充这些线时，不想将线索分开。但是，您的代码对此有所帮助。这是我所拥有的最终代码，可以使您适应每个唯一的Linestring角度。

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"])

Answer 3

您可能有兴趣使用空间网络的 sfnetworks 软件包。它包括函数 edge_azimuth（）来计算网络中边缘的方位角。

请注意，OSM数据可能需要进行一些预处理，然后才能形成一个干净的可路由网络。参见有关预处理和清洁的sfnetworks小插图：

preprex：

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）