加拿大道路网Shapefile导入R遇内存不足及几何列缺失问题求助

一、先修复分块读取的报错问题 #

你用st_read时添加的layer_options参数干扰了sf对几何列的识别，直接去掉多余参数，用基础分块语法即可：

library(sf)
# 读取前100条要素
test <- st_read(dsn = temp_dir, layer = "lrnf000r22a_e", n_max = 100)
# 跳过前100条，读取接下来的100条
test_skip <- st_read(dsn = temp_dir, layer = "lrnf000r22a_e", skip = 100, n_max = 100)

GEOMETRY=AS_WKT会强制将几何转成文本格式，反而导致sf无法识别几何列；FEATURE_TYPE也无需手动指定，sf会自动识别Shapefile的几何类型。

二、解决内存不足的核心方案 #

224万条要素的路网文件远超过你1.46GiB的可用内存，直接全量导入不现实，以下是几个实用方案：

1. 只导入需要的区域/字段 #

• 空间筛选：如果你只需要某省/城市的道路，先获取该区域的边界（比如加拿大省界Shapefile），用GDAL的空间筛选参数提前裁剪路网：

# 假设你已有安大略省边界的sf对象ontario_boundary
boundary_wkt <- st_as_text(st_union(ontario_boundary))
# 只导入安省范围内的道路
ontario_roads <- st_read(
  dsn = temp_dir, 
  layer = "lrnf000r22a_e",
  layer_options = paste0("SPATIAL_FILTER=", boundary_wkt)
)

• 字段筛选：路网的21个字段中多数可能无用，用query参数指定只读取必要字段：

# 仅读取几何列和道路类型字段（假设字段名为ROAD_TYPE）
roads_light <- st_read(
  dsn = temp_dir, 
  layer = "lrnf000r22a_e",
  query = "SELECT geometry, ROAD_TYPE FROM lrnf000r22a_e"
)

2. 用GDAL虚拟矢量（VRT）分块处理 #

创建虚拟矢量文件拆分大Shapefile，无需一次性加载全量数据：

library(sf)
# 按要素数量分块创建VRT文件
gdal_utils(
  util = "buildvrt",
  source = file.path(temp_dir, "lrnf000r22a_e.shp"),
  destination = file.path(temp_dir, "roads.vrt"),
  options = c("-featurelimit", "100000") # 每10万条要素为一个块
)
# 分块读取VRT
chunk1 <- st_read(file.path(temp_dir, "roads.vrt"), n_max = 100000)
chunk2 <- st_read(file.path(temp_dir, "roads.vrt"), skip = 100000, n_max = 100000)

3. 直接用在线API计算驾驶距离（推荐） #

如果你的最终目的是计算两点间驾驶距离，完全无需导入路网。用osrm包调用开放路线服务API（需申请免费密钥）：

library(osrm)
# 配置OpenRouteService端点（密钥需自行官网申请）
options(osrm.server = "https://api.openrouteservice.org/v2/")
options(osrm.profile = "driving-car")
# 定义起点（多伦多）和终点（蒙特利尔）坐标
origin <- st_sfc(st_point(c(-79.38, 43.65)), crs = 4326)
destination <- st_sfc(st_point(-73.57, 45.50), crs = 4326)
# 计算驾驶路线与距离
distance_result <- osrmRoute(origin, destination, returnclass = "sf")
# 提取距离（单位：米）
distance_m <- st_length(distance_result)

4. 优化R内存使用 #

• 关闭其他占用内存的程序，释放更多空间；
• 手动清理内存垃圾：

gc() # 运行后查看可用内存变化

• 调高R的内存限制（根据总内存调整）：

memory.limit(size = 6000) # 单位：MB，设置为6GiB左右

三、总结 #

优先选择在线API计算距离的方案，完全避开大文件导入问题；若必须导入路网，先通过区域/字段筛选减少数据量，再配合分块读取。分块时不要添加干扰几何识别的参数。

内容的提问来源于stack exchange，提问作者stats_noob