r - 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 数据源。

没有查询:

> 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 选项仅在类 DBIObject 的文档中提及，对于“默认 S3”方法，没有 query 参数，因此您的查询字符串被传递到 ... 参数，该参数最终被传递给 st_as_sf ，然后在稍后的某个时间点它会在工作中抛出 Spanner 。

可能有一种解决方法，但一种解决方案是创建一个包含 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 几何列是什么...