st_read() 中的查询问题

Question

我是 sf 包的新手，并尝试读取 shapefile 并根据查询对其进行子集化。在这里，我使用了 sf_read()

  load <- st_read(dsn = "~Data", layer = "CBSA_MetroDiv", 
            query = "select * from CBSA_MetroDiv limit 3;")

但它抛出一个错误

Reading layer `CBSA_MetroDiv' from data source `\Data' using driver `ESRI Shapefile'

Error in st_sf(x, ..., agr = agr, sf_column_name = sf_column_name) : 
no simple features geometry column present

有人可以指导我解决这个问题吗？

Answer 1

更新：query现已实施

该query选项现在应该适用于 GDAL/OGR 数据源。

没有query：

> s = st_read(f)
Reading layer `uk_LAD_may_2020_with_insets_v1.1' from data source 
  `/nobackup/rowlings/Downloads/uk_LAD_may_2020_with_insets_v1.1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 464 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87544.11 ymin: 5344.479 xmax: 980803.6 ymax: 1220302
Projected CRS: OSGB 1936 / British National Grid

获得 464 个特征。通过查询（我引用了图层名称，因为它里面有点，您可能不需要）：

> s = st_read(f, query="select * from \"uk_LAD_may_2020_with_insets_v1.1\" where scale = 3")
Reading query `select * from "uk_LAD_may_2020_with_insets_v1.1" where scale = 3' from data source `/nobackup/rowlings/Downloads/uk_LAD_may_2020_with_insets_v1.1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 52 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87544.11 ymin: 146818.5 xmax: 949029.8 ymax: 986035
Projected CRS: OSGB 1936 / British National Grid

只有 52 个功能。

以前的答案：sf仅适用于较旧的软件包版本：

该query选项仅在 class 的文档中提到DBIObject，对于“默认 S3”方法，没有query参数，因此您的查询字符串将被传递到...最终传递给的参数中st_as_sf，然后在稍后的某个时刻它会在作品。

可能有一种解决方法，但一种解决方案是创建一个包含 SQL 的虚拟数据集文件。例如，我有一个法国邮政区域的 shapefile，这是一个名为的虚拟数据集文件，filter.vrt该文件应用 SQL 选择：

<OGRVRTDataSource>
    <OGRVRTLayer name="points">
        <SrcDataSource relativeToVRT="1">codes_postaux_region.shp</SrcDataSource>
        <SrcSQL>select * from codes_postaux_region where POP2010 > 20000</SrcSQL>
 </OGRVRTLayer>
</OGRVRTDataSource>

使用纯文本编辑器为 shapefile 和 SQL 创建一个类似的文件，然后读取它。在这里你可以看到，如果我读取 shapefile，我会得到 6048 个特征，但当我读取虚拟数据文件时，我只会得到 707 个特征：

> fr = st_read("./codes_postaux_region.shp",quiet=TRUE)
> nrow(fr)
[1] 6048
> fr = st_read("./filter.vrt",quiet=TRUE)
> nrow(fr)
[1] 707

您可能需要在已过滤的数据集上设置坐标系，可以在读入该坐标系（如果您知道）之后，也可以通过另一个 VRT 文件参数来设置。

可能值得 ping Edzer，看看是否st_read可以实现 shapefile 的 SQL，或者我是否遗漏了某些内容。我觉得有一种方法可以告诉st_read几何列是什么......