将线 shapefile 转换为栅格图层并保留值答案

【问题标题】：Convert a line shapefile to a raster layer and preserve values将线 shapefile 转换为栅格图层并保留值
【发布时间】：2021-05-07 11:51:17
【问题描述】：

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

（代码取自this blog post）

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 转换为保留道路类型的栅格文件，以便计算成本路径。

【问题讨论】：

查看光栅包rdocumentation.org/packages/raster/versions/3.4-5/topics/…中的光栅化函数。或者，您也可以使用更快的版本 - rdocumentation.org/packages/fasterize/versions/1.0.3/topics/…
您可能还想检查sfnetworks 中的st_network_cost() 函数，该函数与您提到的博客文章链接。它不提供光栅文件，但也许它仍然对您有用。
谢谢两位——这可能超出了本文的范围，但我的主要目标是在考虑道路质量、拓扑等因素的同时找出两点之间的最佳路径。使用 shapefile 和st_network_cost 建议的功能，我无法考虑非公路路线，所以我认为这是一个基于栅格的问题。现在的问题变成了如何将道路、拓扑等纳入这个成本函数。
@PabloHerrerosCantis 加速功能不适用于线条，仅适用于多边形。请参阅我的答案，以使用 stars 包更快地栅格化线条

标签： r raster spatial sf

【解决方案1】：

您可能已经注意到，对于大型数据集，光栅化非常慢。我创建了一个环绕 stars::st_rasterize 的函数，以便它接受一个 RasterLayer 和 sf 并返回一个 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

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

【讨论】：

【解决方案2】：

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)

【讨论】：

【解决方案3】：

转换为空间线对象：

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 命令来计算成本路径。

【讨论】：