【问题标题】：Issue with query in st_read()st_read() 中的查询问题
【发布时间】：2019-02-04 19:11:46
【问题描述】：

我是 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

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

【问题讨论】：

在没有其他所有参数的情况下是否可以工作，例如load <- st_read("path/to/your/shape/file.shp")
是的。它适用于查询，但在我添加查询时抛出错误。我的目标是只加载形状文件的子集，而不是将整个形状文件加载到内存中。
已经做了一些测试，在这里问：github.com/r-spatial/sf/issues/834 否则检查答案中的VRT解决方法。

标签： r sp sf

【解决方案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 选项仅在 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 几何列是什么...

【讨论】：

更新：query 现已实施

上一个答案：仅适用于较旧的 sf 软件包版本：

更新：`query` 现已实施

上一个答案：仅适用于较旧的 `sf` 软件包版本：