将线形状文件转换为栅格图层并保留值

Question

我有一个 sf 对象，其中包含道路几何和道路类型（例如二级公路、高速公路、自行车道等）：

library(osmdata)
library(tidyverse)

muenster <- 
  opq(bbox =  c(7.61, 51.954, 7.636, 51.968)) %>% 
  add_osm_feature(key = 'highway') %>% 
  osmdata_sf() %>% 
  osm_poly2line()

muenster_center <- muenster$osm_lines %>% 
  dplyr::select(highway)

我想将此 shapefile 转换为保留道路类型的光栅文件，以便我可以计算成本路径。

Answer 1

转换为空间线对象：

muenster_center <- as(muenster_center, "Spatial")

将高速公路变量转化为因子：

muenster_center@data$highway <- as.factor(muenster_center@data$highway)

创建光栅（ncol + nrow定义光栅的粒度）

r <- raster(muenster_center, ncol=1000, nrow=1000)    
m_ras <- rasterize(muenster_center, r, "highway")

运行 plot(m_ras) 生成下图：

然后您可以使用 gdistance 包中的 shortestPath 命令来计算成本路径。

Answer 2

您可能已经注意到，对于大型数据集，栅格化速度非常慢。我创建了一个环绕 stars::st_rasterize 的函数，因此它需要一个 RasterLayer 和 sf 以及 returns 一个 RasterLayer

rasterizeLine <- function(sfLine, rast, value){
  # rasterize roads to template
  tmplt <- stars::st_as_stars(sf::st_bbox(rast), nx = raster::ncol(rast),
                              ny = raster::nrow(rast), values = value)
  
  rastLine <- stars::st_rasterize(sfLine,
                                  template = tmplt,
                                  options = "ALL_TOUCHED=TRUE") %>%
    as("Raster")
  
  return(rastLine)
}

library(osmdata)
library(tidyverse)
library(sf)

muenster <- 
  opq(bbox =  c(7.61, 51.954, 7.636, 51.968)) %>% 
  add_osm_feature(key = 'highway') %>% 
  osmdata_sf() %>% 
  osm_poly2line()

muenster_center <- muenster$osm_lines %>% 
  select(highway) %>% 
  # Note you have to turn it into a factor for stars to use that column and it
  # has to be the first column
  mutate(highway = as.factor(highway))

r <- raster::raster(muenster_center, ncol=1000, nrow=1000)

system.time({
  muenster_rast <- rasterizeLine(muenster_center, r, NA)
})
#   user  system elapsed 
#   0.09    0.07    0.15 

system.time({
  muenster_rast2 <- raster::rasterize(muenster_center, r, "highway")
})
#   user  system elapsed 
# 207.56    0.06  208.75

raster::freq(muenster_rast - muenster_rast2)
#      value  count
# [1,]    -2      1
# [2,]     0  63676
# [3,]    NA 936323

两者的输出看起来非常相似，但使用星标的版本要快得多！

Answer 3

terra 是这样的

library(osmdata)

muenster <- 
  opq(bbox =  c(7.61, 51.954, 7.636, 51.968)) |>
  add_osm_feature(key = 'highway') |>
  osmdata_sf() |> 
  osm_poly2line()

muenster_center <- muenster$osm_lines["highway"]

library(terra)

v <- vect(muenster_center)
r <- rast(v, resolution=1/12000)

x <- rasterize(v, r, "highway")
plot(x)

如果你需要一个 RasterLayer 你可以做

r <- raster::raster(x)

将线形状文件转换为栅格图层并保留值

Convert a line shapefile to a raster layer and preserve values

r

raster

spatial

sf