ArcGIS REST API query by distance from point

Question

I am using the arcpullr package to query GIS data hosted on an ArcGIS server. I am able to download a subset of the data using attributes in my query, but I am not able to get the spatial query get_layer_by_point() working.

Ultimately what I am hoping to do is extract the coordinates of the nearest point on the line for use in another process. I'm open to using another package or building the query outside of a package.

Here is an example. I am able to extract the line using its Permanent_Identifier. I manually looked up this value, but in future cases would prefer to query the line using its proximity to a point.

In this example, I know the point and the line both exist and are close to each other.

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')

^{Created on 2021-09-24 by the reprex package (v2.0.0)}

But I am not able to do so by spatial query relative to the point. The API documentation is here (note the distance parameter) and code I have tried is below

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

score 1 · Accepted Answer · answered Sep 27 '21 at 18:45

Worked out a solution by leaning on the sf package rather than arcpullr as @Jindra Lacko suggested. If I buffer the point using st_buffer I am able to query using the buffer polygon inside of get_layer_by_poly().

To snap the point to the line I followed the answer provided here: https://stackoverflow.com/a/51300037/12400385

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 https://stackoverflow.com/a/51300037/12400385
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)

^{Created on 2021-09-27 by the reprex package (v2.0.0)}

score 0 · Answer 2 · answered Sep 24 '21 at 19:59

0

If you can get the point and line geometry working in R you should be able to produce the desired result using sf::st_nearest_points(). Your code is not exactly reproducible, but have a look at the documentation, it should get you started.

The function supports point / line geometry pairs, and it will produce a linestring of two points: one should be your original point (disregard it) and the other the closest point on your line.

You may find it helpful to sf::st_cast() the output to points type geometry, depending on your application.

answered Sep 24 '21 at 19:59

Jindra Lacko

7,814
3
22
44

Thanks for the tip, `st_nearest_point` will definitely be useful. But my issue is I don't want to have to manually look up the `Permanent_Identifier`. So I would like to query the line by proximity to the point. – nniloc Sep 24 '21 at 20:18
Oki, I must have misunderstood. What about `sf::st_nearest_feature()` then? Can you get a data frame of all possible features to choose the nearest one? – Jindra Lacko Sep 25 '21 at 08:17
I am hoping to avoid downloading the entire layer first, but perhaps downloading once and caching is not such a bad solution. I worked out a solution using `sf` as you suggested, and I found a few hints on snapping a point to a line on SO and GH issues. I appreciate the feedback you gave here – nniloc Sep 27 '21 at 18:48

ArcGIS REST API query by distance from point

2 Answers2