0

我曾使用 OSM 数据库中的中心性计算。我按照说明（链接）。

但是，我想导出带有中心值的街道网络的 shapefile。我是 R 新手。所以我不确定哪个部分用作要导出的对象。

我尝试了2种方法

  1. 写OGR

     writeOGR("edges", dsn="210217", layer="edges", driver = "ESRI Shapefile")

    那么它总是显示 SRI Shapefile") : inherits(obj, "Spatial") is not TRUE

  2. st_write

     st_write(graph, "bkk.shp")
 sf::st_write(graph, "bkk.shp", driver = "ESRI Shapefile")
 Error in UseMethod("st_write") :
 no applicable method for 'st_write' applied to an object of class 
 "c('tbl_graph', 'igraph')"

-------------------------------- 代码我将从这部分开始----

我想导出这张地图的 shapefile：

在此处输入图像描述

writeVECT(
  SDF = bangkok_center, 
  vname = 'bangkok_center', 
  v.in.ogr_flags = 'overwrite'
)

execGRASS("g.proj", flags = c("c", "quiet"), proj4 = proj4)
execGRASS(
  cmd = 'v.clean', 
  input = 'bangkok_center', 
  output = 'bangkok_cleaned',        
  tool = 'break', 
  flags = c('overwrite', 'c')
)

edges <-bangkok_center %>%
  mutate(edgeID = c(1:n()))
edges

nodes <- edges %>%
  st_coordinates() %>%
  as_tibble() %>%
  rename(edgeID = L1) %>%
  group_by(edgeID) %>%
  slice(c(1, n())) %>%
  ungroup() %>%
  mutate(start_end = rep(c('start', 'end'), times = n()/2))

nodes

nodes <- nodes %>%
  mutate(xy = paste(.$X, .$Y)) %>% 
  mutate(nodeID = group_indices(., factor(xy, levels = unique(xy)))) %>%
  select(-xy)
nodes

source_nodes <- nodes %>%
  filter(start_end == 'start') %>%
  pull(nodeID)

target_nodes <- nodes %>%
  filter(start_end == 'end') %>%
  pull(nodeID)

edges = edges %>%
  mutate(from = source_nodes, to = target_nodes)

edges

nodes <- nodes %>%
  distinct(nodeID, .keep_all = TRUE) %>%
  select(-c(edgeID, start_end)) %>%
  st_as_sf(coords = c('X', 'Y')) %>%
  st_set_crs(st_crs(edges))

nodes

graph = tbl_graph(nodes = nodes, edges = as_tibble(edges), directed = FALSE)
graph

sf_to_tidygraph = function(x, directed = TRUE) {
  
  edges <- x %>%
    mutate(edgeID = c(1:n()))
  
  nodes <- edges %>%
    st_coordinates() %>%
    as_tibble() %>%
    rename(edgeID = L1) %>%
    group_by(edgeID) %>%
    slice(c(1, n())) %>%
    ungroup() %>%
    mutate(start_end = rep(c('start', 'end'), times = n()/2)) %>%
    mutate(xy = paste(.$X, .$Y)) %>% 
    mutate(nodeID = group_indices(., factor(xy, levels = unique(xy)))) %>%
    select(-xy)
  
  source_nodes <- nodes %>%
    filter(start_end == 'start') %>%
    pull(nodeID)
  target_nodes <- nodes %>%
    filter(start_end == 'end') %>%
    pull(nodeID)
  edges = edges %>%
    mutate(from = source_nodes, to = target_nodes)
  
  nodes <- nodes %>%
    distinct(nodeID, .keep_all = TRUE) %>%
    select(-c(edgeID, start_end)) %>%
    st_as_sf(coords = c('X', 'Y')) %>%
    st_set_crs(st_crs(edges))
  
  tbl_graph(nodes = nodes, edges = as_tibble(edges), directed = directed)
  
}
sf_to_tidygraph(bangkok_center, directed = FALSE)

graph <- graph %>%
  activate(edges) %>%
  mutate(length = st_length(geometry))
graph

graph %>%
  activate(edges) %>%
  as_tibble() %>%
  st_as_sf() %>%
  group_by(highway) %>%
  summarise(length = sum(length))



ggplot() +
  geom_sf(data = graph %>% activate(edges) %>% as_tibble() %>% st_as_sf()) + 
  geom_sf(data = graph %>% activate(nodes) %>% as_tibble() %>% st_as_sf(), size = 0.5)

graph <- graph %>%
  activate(nodes) %>%
  mutate(degree = centrality_degree()) %>%
  mutate(betweenness = centrality_betweenness(weights = length)) %>%
  activate(edges) %>%
  mutate(betweenness = centrality_edge_betweenness(weights = length))

graph

st_write(graph, "bkk.shp")
sf::st_write(graph, "bkk.shp", driver = "ESRI Shapefile")

ggplot() +
  geom_sf(data = graph %>% activate(edges) %>% as_tibble() %>% st_as_sf(), aes(col = betweenness, size = betweenness)) +
  scale_colour_viridis_c(option = 'inferno') +
  scale_size_continuous(range = c(0,4))

ggplot() +
  geom_sf(data = graph %>% activate(edges) %>% as_tibble() %>% st_as_sf(), col = 'grey50') + 
  geom_sf(data = graph %>% activate(nodes) %>% as_tibble() %>% st_as_sf(), aes(col = degree, size = degree)) +
  scale_colour_viridis_c(option = 'inferno') +
  scale_size_continuous(range = c(0,4))
4

1 回答 1

0

错误似乎对象graph不是空间对象，而是由两组不同观察值组成的网络类型图 1. 节点或点；2.边缘或链接。

尝试先将graph的边和节点转换为sf格式，然后将它们导出到 shapefile。IE：

sf_graph_nodes <- graph %>% activate(nodes) %>% as_tibble() %>% st_as_sf()
sf_graph_edges <- graph %>% activate(edges) %>% as_tibble() %>% st_as_sf()

st_write(sf_graph_nodes, "nodes_bkk.shp")
st_write(sf_graph_edges, "edges_bkk.shp")
于 2021-02-18T23:40:21.253 回答