ArcGIS REST API 按距点的距离查询

Question

我正在使用 arcpullr package 查询托管在 ArcGIS 服务器上的 GIS 数据。我可以使用查询中的属性下载数据的子集，但我无法使空间查询 get_layer_by_point() 正常工作。

最终我希望做的是提取直线上最近点的坐标以用于另一个过程。我愿意使用其他包或在包外构建查询。

这是一个例子。我能够使用其 Permanent_Identifier 提取该行。我手动查找了这个值，但在未来的情况下，我更愿意使用它与某个点的接近度来查询该线。

在这个例子中，我知道点和线都存在并且彼此靠近。

library(arcpullr)
#> Loading required package: sf
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
  
pt_url <- "https://geodataservices.wdfw.wa.gov/arcgis/rest/services/ApplicationServices/FP_Sites/MapServer/1"
pt_sf <- get_spatial_layer(pt_url, where = "SiteId = '991057'")
#> Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
#> deprecated. It might return a CRS with a non-EPSG compliant axis order.
  

# works to extract this reach using its Identifier code
line_url <- "https://fortress.wa.gov/ecy/gisprod/arcgis/rest/services/NHD/NHD_Cache/MapServer/1"
line_sf <- get_spatial_layer(line_url, where = "Permanent_Identifier  = '165793370'")


plot(st_geometry(line_sf))
plot(st_geometry(pt_sf), add = TRUE, col = 'red')

^{由 reprex package (v2.0.0)}

于 2021-09-24 创建

但是我无法通过相对于该点的空间查询来做到这一点。 API documentation 在这里（注意 distance 参数），我试过的代码在下面

line_sf2 <- get_layer_by_point(line_url, pt_sf, where = sql_where(distance=500, units="feet"))
#> Warning in get_esri_features(query_url, out_fields, where, token, ...): No
#> records match the search critera

Answer 1

如果您可以在 R 中使用点和线几何图形，您应该能够使用 sf::st_nearest_points() 生成所需的结果。您的代码不能完全重现，但请查看文档，它应该可以帮助您入门。

该函数支持点/线几何对，它会生成一个包含两个点的线串：一个应该是您的原始点（忽略它），另一个是您线上最近的点。

根据您的应用程序，您可能会发现 sf::st_cast() 输出到点类型几何图形很有帮助。

Answer 2

依靠 sf 包而不是 @Jindra Lacko 建议的 arcpullr 来制定解决方案。如果我使用 st_buffer 缓冲点，我可以使用 get_layer_by_poly().

内的缓冲多边形进行查询

为了将点捕捉到我遵循此处提供的答案的线：

library(arcpullr)
#> Loading required package: sf
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1


pt_url <- "https://geodataservices.wdfw.wa.gov/arcgis/rest/services/ApplicationServices/FP_Sites/MapServer/1"
pt_sf <- get_spatial_layer(pt_url, where = "SiteId = '991057'")
#> Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
#> deprecated. It might return a CRS with a non-EPSG compliant axis order.


# works to extract this reach using its Identifier code
line_url <- "https://fortress.wa.gov/ecy/gisprod/arcgis/rest/services/NHD/NHD_Cache/MapServer/1"


pt_buff <- st_buffer(pt_sf, units::set_units(500, 'feet'))
streams <- get_layer_by_poly(line_url, pt_buff, sp_rel = "esriSpatialRelIntersects")


# copied from 
st_snap_points = function(x, y, max_dist = 1000) {
  
  if (inherits(x, "sf")) n = nrow(x)
  if (inherits(x, "sfc")) n = length(x)
  
  out = do.call(c,
                lapply(seq(n), function(i) {
                  nrst = st_nearest_points(st_geometry(x)[i], y)
                  nrst_len = st_length(nrst)
                  nrst_mn = which.min(nrst_len)
                  if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
                  return(st_cast(nrst[nrst_mn], "POINT")[2])
                })
  )
  return(out)
}


pt_snap <- st_snap_points(pt_sf, streams)

plot(st_geometry(pt_buff))
plot(streams["Permanent_Identifier"], add = TRUE)
plot(st_geometry(pt_sf), add = TRUE)
plot(st_geometry(pt_snap), col = 'red', add = TRUE)

^{由 reprex package (v2.0.0)}

于 2021-09-27 创建

ArcGIS REST API 按距点的距离查询

ArcGIS REST API query by distance from point

r

gis

sf

arcpullr