使用{sf}包查询与多边形相交的GeoPackage要素报错求助

错误根源 #

1. 空间函数调用不符合GeoPackage规范：你使用的st_polygon调用方式不是GeoPackage SQL支持的格式。GeoPackage遵循OGC空间函数规范，构造多边形需要用ST_GeomFromText（传入标准WKT）或ST_MakeEnvelope（专门用于边界框）。
2. WKT格式错误：你传入的坐标串没有包裹在正确的POLYGON WKT结构中。标准POLYGON WKT必须是POLYGON((x1 y1, x2 y2, ..., x1 y1))的格式，多层括号是必须的。

修复后的代码 #

推荐方案：使用ST_MakeEnvelope（更高效简洁） #

ST_MakeEnvelope直接接收边界框的四个坐标值和坐标系EPSG代码，无需手动拼接WKT，是处理边界框相交查询的最优选择：

box <- st_read("file1.gpkg", quiet=T) %>% st_bbox()
layer_name <- st_layers("file2.gpkg")$name

# 获取目标图层的坐标系EPSG代码（避免硬编码出错）
target_crs <- st_crs(st_read("file2.gpkg", quiet=T))$epsg

my_query <- glue("SELECT * FROM {layer_name} WHERE ST_Intersects(geom, ST_MakeEnvelope({box$xmin}, {box$ymin}, {box$xmax}, {box$ymax}, {target_crs}))")

st_read("file2.gpkg", query = my_query)

备选方案：使用ST_GeomFromText构造标准WKT #

如果需要手动构造多边形WKT，确保格式符合标准：

box <- st_read("file1.gpkg", quiet=T) %>% st_bbox()
layer_name <- st_layers("file2.gpkg")$name

target_crs <- st_crs(st_read("file2.gpkg", quiet=T))$epsg

# 生成符合规范的POLYGON WKT
polygon_wkt <- glue("POLYGON(({box$xmin} {box$ymin}, {box$xmax} {box$ymin}, {box$xmax} {box$ymax}, {box$xmin} {box$ymax}, {box$xmin} {box$ymin}))")
my_query <- glue("SELECT * FROM {layer_name} WHERE ST_Intersects(geom, ST_GeomFromText('{polygon_wkt}', {target_crs}))")

st_read("file2.gpkg", query = my_query)

关键注意事项 #

• 必须确保查询中使用的坐标系EPSG代码与目标GeoPackage图层的坐标系完全一致，否则空间相交判断会出错或直接报错。
• GeoPackage的SQL空间函数名称通常是大写的（如ST_Intersects），虽然部分引擎支持小写，但统一使用大写能避免兼容性问题。

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