---
title: "Example use cases to get data from web feature or web coverage services."
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{wfs_wcs}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(sf)
library(dplyr)
library(inbospatial)
```

# Introduction

A WFS or Web Feature Service allows you to download vector spatial data from the web.
A WCS or Web Coverage Service is similar, but is used for raster spatial data.
In this vignette, we will show how to the functions `get_feature_wfs()` and `get_coverage_wcs()` can be used to get data into your R session from these services.

For a more in-depth explanation about WFS services, we refer to [the INBO tutorial about WFS](https://tutorials.inbo.be/tutorials/spatial_wfs_services/){target="_BLANK"}.
The tutorial takes a more hands-on coding approach to getting data from a WFS, whereas in this vignette much of the coding stuff is hidden away in the `get_feature_wfs()` function.

# WFS examples

## Simple use case without filter

As a test case, we will query municipal boundaries of Flanders from the "Provisional Reference Municipal Boundaries" (`VRBG`, `Voorlopig ReferentieBestand Gemeentegrenzen`), which is available via an API on the `vlaanderen.be` website.

```{r provinces}
wfs_vrbg <- "https://geo.api.vlaanderen.be/VRBG/wfs"

provinces <- get_feature_wfs(
  wfs = wfs_vrbg,
  layername = "VRBG:Refprv",
  crs = "EPSG:31370"
)

provinces
```

You can achieve (almost) the same result using a combination of `sf::read_sf()` and `sf::st_transform()`.

```{r}
provinces2 <- read_sf(paste0("WFS:", wfs_vrbg),
  layer = "VRBG:Refprv"
)

waldo::compare(provinces, provinces2,
  x_arg = "inbospatial", y_arg = "sf",
  max_diffs = 8
)
```

Both methods yield the same result.
Note however that for other examples, you may notice substantial differences, most notably for field specifications.
The reason for such differences is that `get_feature_wfs()` retrieves whatever information that is identified from the WFS request - in the same way as you would paste the request in the browser.
This result is then saved temporarily to a `.GML` file and read with `sf::read_sf()`.
Passing the WFS directly to the `dsn` argument of `sf::read_sf`, on the other hand, will translate the request to a form that will pass through [the GDAL WFS driver](https://gdal.org/en/latest/drivers/vector/wfs.html), which handles field specifications by reading them from a `DescribeFeatureType` request:

> At the first opening, the content of the result of the `GetCapabilities` request will be appended to the file, so that it can be cached for later openings of the dataset. The same applies for the `DescribeFeatureType` request issued to discover the field definition of each layer.

## Apply a CQL-Filter to filter on a field

CQL, the "Common Query Language" ([see here, for example](https://docs.geoserver.org/stable/en/user/tutorials/cql/cql_tutorial.html)), provides a way to filter WFS data.
The syntax of [CQL is comparable to SQL](https://www.webcomand.com/docs/language/cql/cql_vs_sql/), yet there are differences in data types and structures.

You can provide and apply filters in `get_feature_wfs()` via the `cql_filter` argument.

```{r prov}
wfs_vrbg <- "https://geo.api.vlaanderen.be/VRBG/wfs"

prov <- get_feature_wfs(
  wfs = wfs_vrbg,
  layername = "VRBG:Refprv",
  cql_filter = "NAAM='West-Vlaanderen'"
)

prov
```

Here, we could have also used `sf::read_sf()` in combination with an [OGR SQL query](https://gdal.org/en/latest/user/ogr_sql_dialect.html):

```{r}
prov2 <- read_sf(paste0("WFS:", wfs_vrbg),
  query = "SELECT * FROM Refprv WHERE NAAM='West-Vlaanderen'"
)

waldo::compare(prov, prov2)
```

Both methods yield the same result.

## Request a subset of fields

Instead of querying the complete data table, one can select a subset of variables (fields).
The function `get_feature_wfs()` provides the `property_name` keyword for it, which accepts a string of comma-separated field names.
Note that field names are capitalized on this particular geoserver.

```{r prov-fields}
wfs_vrbg <- "https://geo.api.vlaanderen.be/VRBG/wfs"

prov_fields <- get_feature_wfs(
  wfs = wfs_vrbg,
  layername = "VRBG:Refprv",
  crs = "EPSG:31370",
  property_name = paste("NAAM", "NUTS2", "SHAPE", sep = ",")
)

prov_fields
```

Again, we can use `sf::read_sf()` in combination with an [OGR SQL query](https://gdal.org/en/latest/user/ogr_sql_dialect.html):

```{r}
prov_fields_sf <- read_sf(paste0("WFS:", wfs_vrbg),
  query = "SELECT NAAM, NUTS2 FROM Refprv"
)

waldo::compare(prov_fields, prov_fields_sf)
```

## Meteorological data example

Another example of a WFS API is the Belgian Royal Meteorological Institute (KMI, `Koninklijk Meteorologisch Instituut van België`).
Details can be [found here](https://opendata.meteo.be/service/wfs?request=GetCapabilities), [and here](https://publish.geo.be/geonetwork/f0ow2say/api/records/RMI_AWS_WFS).

```{r stations}
kmi <- "https://opendata.meteo.be/service/wfs"

kmi_stations <- get_feature_wfs(
  wfs = kmi,
  layer = "synop:synop_station"
)

kmi_stations
```

This is a list of all KMI stations.
You can plot it with `leaflet` or `mapview` to see the stations on a map and learn that -for instance- Ukkel has code 6447.
So let's extract the maximum daily temperatures for Ukkel starting from the year 2000.

```{r kmi, error=TRUE}
kmi_synop <- get_feature_wfs(
  wfs = kmi,
  layername = "synop:synop_data",
  property_name = paste(
    "code", "timestamp", "temp_max", "the_geom",
    sep = ","
  ),
  cql_filter = "((code = 6447) AND (temp_max IS NOT NULL) AND
    (timestamp >= '2000-01-01 00:00:00.000'))"
)

kmi_synop
```


## Spatial filtering using a bounding box

Few analyses will require a whole countries' data, even if the country is as small as Belgium.
Hence, it might be useful by restricting the data you query with a bounding box.
This works as follows, using a `bbox` with the elements `xmin`, `ymin`, `xmax`, and `ymax`.

```{r bwk-bbox}
wfs_bwk <-
  "https://geo.api.vlaanderen.be/BWK/wfs"

bwk_bbox <- get_feature_wfs(
  wfs = wfs_bwk,
  layername = "BWK:Bwkhab",
  crs = "EPSG:31370",
  bbox = c(xmin = 142600, ymin = 153800, xmax = 146000, ymax = 156900)
)

bwk_bbox
```

This returns a `MULTISURFACE` geometry type.
We can use a couple of `sf` functions to cast this to regular polygons:

```{r}
bwk_bbox %>%
  st_cast("GEOMETRYCOLLECTION") %>%
  st_collection_extract("POLYGON") %>%
  st_geometry() %>%
  plot()
```

This can also be done using the `wkt_filter` argument from `sf::read_sf`.

```{r eval=FALSE}
bwk_bbox_sf <- read_sf(
  paste0("WFS:", wfs_bwk),
  layer = "BWK:Bwkhab",
  wkt_filter = st_as_text(
    st_as_sfc(
      st_bbox(
        c(xmin = 142600, ymin = 153800, xmax = 146000, ymax = 156900),
        crs = st_crs(31370)
      )
    )
  )
)
```


## Spatial filters using spatial binary predicates

In addition to the filters above, you can use logical operations to restrict your area of interest.
For example, features can be selected by intersection with another geographical object.

### Intersects a point

As a proof of concept, we can directly query soil parameters at a specific location.

```{r}
wfs_soiltypes <-
  "https://www.dov.vlaanderen.be/geoserver/bodemkaart/bodemtypes/wfs"

x_lam <- 173995
y_lam <- 212093

soil_type <- get_feature_wfs(
  wfs = wfs_soiltypes,
  layername = "bodemkaart:bodemtypes",
  crs = "EPSG:31370",
  cql_filter = sprintf(
    "INTERSECTS(geom,POINT(%s %s))",
    x_lam, y_lam
  )
)

soil_type
```

In theory, this should also be possible "the `sf` way".
However, `sf` filters after query (which takes a long time), and the code below does not work correctly.

```{r eval=FALSE}
soil_type_sf <- read_sf(
  paste0("WFS:", wfs_soiltypes),
  query = "SELECT * FROM bodemtypes",
  wkt_filter = sprintf(
    "ST_Intersects(geom,ST_Point(%s, %s))",
    x_lam, y_lam
  )
)

soil_type_sf
```

# WCS examples

Above, we covered Web Feature Services (WFS), which are many kinds of vector data.
Another common geographical data type are grid-like or "raster" data, which you can retrieve via *Web Coverage Services (WCS)*.
Behold, there is a `get_coverage_wcs()` available for some of the web services we regularly use at INBO.

Below the example of a false-colour infrared image.
The three bands of the image are visualized using `terra::plotRGB`.

```{r omz-example}
bbox <- sf::st_bbox(
  c(xmin = 99477, xmax = 99541, ymin = 172580, ymax = 172640),
  crs = sf::st_crs(31370)
)
omz <- get_coverage_wcs(
  wcs = "omz",
  bbox = bbox,
  layername = "OI.OrthoimageCoverage.OMZ.CIR",
  resolution = 1,
  version = "1.0.0"
)
terra::plotRGB(omz)
```

Another example, this time extracting the digital terrain model for the same bounding box:

```{r dtm-example}
dtm <- get_coverage_wcs(
  wcs = "dtm",
  bbox = bbox,
  layername = "EL.GridCoverage.DTM",
  resolution = 1)
terra::plot(dtm)
```

More details on WCS can be found in the `spatial_dhmv_query.Rmd` vignette.