如何使用st_read读取时实现GeoPackage与另一GeoPackage子集相交

你可以利用GeoPackage的空间SQL能力，在st_read的query参数中直接执行空间过滤，只读取和目标边界相交的建筑物数据，完全不需要加载整个大数据集。具体实现步骤如下：

1. 获取目标边界的WKT格式 #

先读取并筛选出你需要的边界，然后转换成WKT（Well-Known Text）格式，用于后续的SQL查询：

library(sf)

# 读取边界数据集并筛选目标区域
target_boundary <- st_read("https://automaticknowledge.org/gb/wards/Manchester_wards.gpkg") %>%
  filter(WD21CD == "E05011363")

# 将目标边界转换为WKT格式（确保坐标系与建筑物数据集一致，这里是EPSG:4326）
target_wkt <- st_as_text(st_geometry(target_boundary))

2. 用空间SQL查询直接读取相交数据 #

构造SQL语句，通过ST_Intersects函数在数据源层面过滤出和目标边界相交的建筑物，再用st_read执行查询：

# 直接读取与目标边界相交的建筑物数据
Buildings_Intersection <- st_read(
  dsn = "https://automaticknowledge.org/gb/buildings/Manchester_buildings.gpkg",
  query = sprintf(
    "SELECT * FROM Manchester_buildings WHERE ST_Intersects(geom, ST_GeomFromText('%s', 4326))",
    target_wkt
  )
)

关键说明 #

• 这里依赖GeoPackage底层的SQLite空间扩展，过滤逻辑在数据源端执行，只会返回符合条件的记录，大幅节省内存和时间。
• 确保两个数据集的坐标系一致，示例中均为WGS84（EPSG:4326），如果不一致，需要在SQL中用ST_Transform转换坐标系，比如：

  ST_Intersects(geom, ST_Transform(ST_GeomFromText('%s', 原EPSG), 目标EPSG))
  

• 若建筑物数据集的几何列不是默认的geom，可以用st_layers("你的GeoPackage路径")查看实际列名并替换。
• 如果目标边界包含多个要素，可以先用ST_Union合并成单个几何对象再执行查询，提升效率。

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