---
title: "Spatial Queries"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  rmarkdown::html_vignette
vignette: > 
  %\VignetteIndexEntry{Spatial Queries}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---

---
title: "Spatial Queries"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to arcpullr}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---
To see more complete package documentation check out:
<a href="https://pfrater.github.io/arcpullr/">
https://pfrater.github.io/arcpullr/</a>
<hr>
```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE
)
library(arcpullr)
library(sf)
sf::sf_use_s2(FALSE)
```

```{r, echo = FALSE}
#<img src='../man/figures/logo.png' width="160" height="180" style="border: none; float: right"/>
```



ArcGIS REST API's may be spatially queried using the get_layer_by\_\* family of 
functions. These functions require a spatial object of 
class `sf` (*i.e.* of the R package [sf:
Simple Features for R](https://r-spatial.github.io/sf/)) and a 
[Spatial Relationship](https://help.arcgis.com/en/webapi/wpf/apiref/ESRI.ArcGIS.Client~ESRI.ArcGIS.Client.Tasks.SpatialRelationship.html)
to be passed to the geometry and `sp_rel` arguments respectively.

The package contains five functions that can be used to perform
spatial queries:

<ul>

<li>

`get_layer_by_line`

<li>

`get_layer_by_point`

<li>

`get_layer_by_polygon`

<li>

`get_layer_by_multipoint`

<li>

`get_layer_by_envelope`

</u >

## URL's for examples 
<details>
<summary
  style="
    display: inline-block;
    padding: 0.5em 1em;
    background-color: #0d6efd;
    color: white;
    border-radius: 0.25rem;
    cursor: pointer;
    font-weight: bold;
    list-style: none;">
<strong>Example Source Data</strong></summary>
```{r url_setup, eval = FALSE}

#WDNR Server
server <- "https://dnrmaps.wi.gov/arcgis/rest/services/"
server2 <- "https://dnrmaps.wi.gov/arcgis2/rest/services/"

#River URL
layer <- "TS_AGOL_STAGING_SERVICES/EN_AGOL_STAGING_SurfaceWater_WTM/MapServer/2"
river_url <- paste0(server2,layer)

#Country URL 
layer <- "DW_Map_Dynamic/EN_Basic_Basemap_WTM_Ext_Dynamic_L16/MapServer/3"
county_url <- paste0(server,layer)

#Trout URL 
layer <- "FM_Trout/FM_TROUT_HAB_SITES_WTM_Ext/MapServer/0"
trout_url <- paste0(server,layer)

#Watershed URL
layer <- "WT_SWDV/WT_Federal_Hydrologic_Units_WTM_Ext/MapServer/0"
watershed_url <- paste0(server,layer)

#get layers for queries
mke_river <- get_spatial_layer(
  river_url, 
  where = "RIVER_SYS_NAME = 'Milwaukee River'"
)

trout_hab_project_pts <- get_spatial_layer(
  trout_url,
  where = "WATERBODYNAMECOMBINED = 'Sugar Creek' and FISCALYEAR = 2017"
)

trout_hab_project_pt <- trout_hab_project_pts[1, ]

# get watershed layer for Cook Creek
cook_creek_ws <- get_spatial_layer(
  watershed_url,
  where = "HUC12_NAME = 'Cook Creek'"
)
```
</details>

\

## get_layer_by_line
The `get_layer_by_line` function uses A LINSESTRING or MULTILINESTRING sf object 
to query an ArcGIS REST API. The below example uses a MULTILINESTRING sf 
object of the Milwaukee River to query the Wisconsin County polygon layer. 

```{r pull_by_line, eval = FALSE, echo = FALSE}
mke_river_counties <- get_layer_by_line(url = county_url, geometry = mke_river)
```

```{r, echo = FALSE}
mke_river_counties <- sf::st_filter(
  wis_counties, 
  mke_river, 
  .pred = sf::st_intersects
)
```

```{r show_by_line, ref.label=c("pull_by_line", "plot_by_line"), eval = FALSE}
```

```{r plot_by_line, echo = FALSE}
plot_layer(mke_river, outline_poly = mke_river_counties)
```

## get_layer_by_point

The `get_layer_by_line` function uses a POINT sf object to query 
an ArcGIS REST API. The below example shows how this can be used to return
which rivers intersect with a trout habitat project on Sugar Creek in southeast
Wisconsin.

```{r pull_by_point, eval = FALSE, echo = FALSE}
trout_stream <- get_layer_by_point(url = river_url, geometry = trout_hab_project_pt)
```

```{r, echo = FALSE, mesage = FALSE}
trout_stream <- sf::st_filter(
  sugar_creek, 
  sf::st_buffer(trout_hab_project_pt, 0.0001), 
  .pred = sf::st_intersects
)
```

```{r show_by_point, ref.label=c("pull_by_point", "plot_by_point"), eval = FALSE}
```

```{r plot_by_point, echo = FALSE}
plot_layer(trout_stream) +
  ggplot2::geom_sf(data = trout_hab_project_pt, color = "red", size = 2)
```




`get_layer_by_point` can also handle multipoint objects. This example shows the
same stream as above with a single point, but now with multiple restoration 
points.

```{r pull_by_multipoint, eval = FALSE, echo = FALSE}
restored_streams <- get_layer_by_point(url = river_url, geometry = trout_hab_project_pts)
```

```{r show_by_multipoint, ref.label=c("pull_by_multipoint", "plot_by_multipoint"), eval = FALSE}
```

```{r, echo = FALSE}
restored_streams <- sugar_creek
```

```{r plot_by_multipoint, fig.height = 7, fig.width = 7, echo = FALSE}
plot_layer(restored_streams) + 
  ggplot2::geom_sf(data = trout_hab_project_pts, color = "blue")
```




## get_layer_by_polygon

The `get_layer_by_line` function uses a POLYGON sf object to query 
an ArcGIS REST API. The below examples shows how this can be used to find
what rivers are within a particular watershed.

```{r, pull_by_poly, eval = FALSE, echo = FALSE}
cook_creek_streams <- `get_layer_by_poly(river_url, cook_creek_ws)
```

```{r show_by_poly, ref.label=c("pull_by_poly", "plot_by_poly"), eval = FALSE}
```

```{r plot_by_poly, echo = FALSE}
plot_layer(cook_creek_streams, cook_creek_ws)
```

## get_layer_by_envelope
The `get_layer_by_envelope` function accepts any sf object to query an
ArcGIS REST API using the sf objects bounding box.  The below example
shows how this is used to query WI's Rivers ArcGIS REST API using a  sf POLYGON 
object of a watershed for a small stream. Note how the results compare
to when this same object is queried using the `get_layer_by_poly` function. 

```{r, pull_by_env, eval = FALSE, echo = FALSE}
cook_creek_env <- get_layer_by_envelope(river_url, cook_creek_ws)

```

```{r show_by_env, ref.label=c("pull_by_env", "plot_by_env"), eval = FALSE}
```

```{r plot_by_env, echo = FALSE}
# example of the envelope to visualize how it spatially queries
example_env <- sf::st_as_sfc(sf::st_bbox(cook_creek_ws))

plot_layer(cook_creek_env, cook_creek_ws) + 
  ggplot2::geom_sf(data = example_env, fill = NA)
```

\
\

## Combining Spatial and SQL Queries
Spatial queries can be combined with SQL statements to further refine queries.

## Spatial Relationship

The `sp_rel` argument can be used to define the spatial relationship
between the two feature classes involved within a spatial query. The
default spatial relationships for the `get_layer_by_poly` function is 
"contains". All other functions default to 
"intersects". 

```{r pull_sp_rel_contains, eval = FALSE, echo = FALSE}
example_poly <- sf_polygon(
  c(-90.62, 43.76),
  c(-90.62, 43.77),
  c(-90.61, 43.77),
  c(-90.61, 43.76),
  c(-90.62, 43.76)
)
poly_streams_contains <- get_layer_by_poly(river_url, example_poly)
```

```{r show_sp_rel_contains, ref.label=c("pull_sp_rel_contains", "plot_sp_rel_contains"), eval = FALSE}
```

```{r plot_sp_rel_contains, echo = FALSE}
plot_layer(poly_streams_contains, outline_poly = example_poly)
```

Using "crosses" returns different records compared to the 
above example (i.e. this returns records when they cross the polygon border).

```{r pull_sp_rel_crosses, eval = FALSE, echo = FALSE}
poly_streams_crosses <- get_layer_by_poly(river_url, example_poly, sp_rel = "crosses")
```

```{r show_sp_rel_crosses, ref.label=c("pull_sp_rel_crosses", "plot_sp_rel_crosses"), eval = FALSE}
```

```{r plot_sp_rel_crosses, echo = FALSE}
plot_layer(poly_streams_crosses, outline_poly = example_poly)
```

<br>

### Lookup Tables
The `sp_rel_lookup` data.frame explains the various types of
spatial relationships available through ArcGIS REST APIs.

```{r, echo = FALSE}
sp_rel_lookup %>%
  DT::datatable()
```
<br>
The `sp_rel_valid` data.frame shows which spatial relationships are
valid with different geometry types being queried and used to do spatial
queries.

```{r, echo = FALSE}
sp_rel_valid %>%
  DT::datatable()

```

### The valid_sp_rel Function 
The `valid_sp_rel`  function can be used to to see which spatial relation
types are applicable to different geometries.

```{r, echo = TRUE}

valid_sp_rel("line","line")

```