diff --git a/NEWS.md b/NEWS.md index cb0bd5e..eba3e0f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,16 @@ # inbospatial 0.0.3 +## Features + +* `get_coverage_wcs()` can now query data from the + Digital Elevation Model, Flanders + (DHMV, ``Digitaal Hoogtemodel Vlaanderen``). + +## Documentation + +* added `spatial_dhmv_query` vignette demonstrating WCS queries + with and beyond `get_coverage_wcs()` + ## Bug fixes * fix parsing of `geoTIFF` part from `mht` files (internal function) diff --git a/R/get_coverage_wcs.R b/R/get_coverage_wcs.R index a28cade..cef0fda 100644 --- a/R/get_coverage_wcs.R +++ b/R/get_coverage_wcs.R @@ -4,7 +4,7 @@ #' from which it is read with `terra::rast()` - if needed reprojected - #' and returned as a `SpatRaster`object #' -#' @param wcs One of `"dtm"`, `"dsm"`, `"omz"`, `"omw"` +#' @param wcs One of `"dtm"`, `"dsm"`, `"omz"`, `"omw"`, `"dhmv"` #' @param bbox An object of class bbox of length 4. #' @param layername Character string; name of the layer #' @param resolution Output resolution in meters @@ -20,7 +20,8 @@ #' - `"omw"`: orthophotomosaic winter images Flanders #' - `"dtm"`: digital terrain model Flanders #' - `"dsm"`: digital surface model Flanders -#' See [metadata Vlaanderen](https://metadata.vlaanderen.be/metadatacenter/srv/dut/catalog.search#/search?keyword=OGC:WCS) for more information # nolint: line_length_linter. +#' - `"dhmv"`: digital elevation model Flanders (contains dtm and dsm data) +#' See [metadata Vlaanderen](https://metadata.vlaanderen.be/srv/eng/catalog.search#/search?any=WCS) for more information # nolint: line_length_linter. #' #' @importFrom sf st_as_sf st_transform st_coordinates #' @importFrom terra rast `res<-` project @@ -46,35 +47,56 @@ #' } #' get_coverage_wcs <- function( - wcs = c("dtm", "dsm", "omz", "omw"), + wcs = c("dtm", "dsm", "omz", "omw", "dhmv"), bbox, layername, resolution, - wcs_crs = "EPSG:4258", + wcs_crs = c("EPSG:4258", "EPSG:31370"), output_crs = "EPSG:31370", bbox_crs = "EPSG:31370", version = c("1.0.0", "2.0.1"), ...) { + # prelim check version <- match.arg(version) + wcs <- tolower(wcs) # case insensitive wcs wcs <- match.arg(wcs) wcs_crs <- match.arg(wcs_crs) bbox_crs <- match.arg(bbox_crs) + + # constrain version | wcs + if (wcs == "dhmv") { + # warn incompatible versions + if ( !( version %in% c("1.0.0", "2.0.1") ) ) { + message("WCS `DHMV` is only compatible with versions `1.0.0` or `2.0.1`. + Consider using `version=\"2.0.1\"`") + } + # recommend crs specification + if (wcs_crs != "EPSG:31370") { + message("WCS `DHMV` only supports CRS Belgian Lambert 72 (`EPSG:31370`). + Consider specifying `wcs_csr=\"EPSG:31370\"`") + } + } + # set url wcs <- switch(wcs, omz = "https://geo.api.vlaanderen.be/oi-omz/wcs", omw = "https://geo.api.vlaanderen.be/oi-omw/wcs", dtm = "https://geo.api.vlaanderen.be/el-dtm/wcs", - dsm = "https://geo.api.vlaanderen.be/el-dsm/wcs" + dsm = "https://geo.api.vlaanderen.be/el-dsm/wcs", + dhmv = "https://geo.api.vlaanderen.be/DHMV/wcs" ) + # data type assertions assert_that(is.character(layername)) assert_that(is.character(output_crs)) assert_that(inherits(bbox, "bbox")) - assert_that(is.numeric(resolution)) + # resolution <=0 will give a `404` + assert_that(is.numeric(resolution) && resolution > 0) + # assemble the bounding box matrix(bbox, ncol = 2, byrow = TRUE) |> as.data.frame() |> st_as_sf(coords = c("V1", "V2"), crs = bbox_crs) |> @@ -83,9 +105,10 @@ get_coverage_wcs <- function( as.vector() -> bbox names(bbox) <- c("xmin", "xmax", "ymin", "ymax") - # build url request + # prepare url request url <- parse_url(wcs) + # variant: version 2.0.1 if (version == "2.0.1") { epsg_code <- str_extract(wcs_crs, "\\d+") url$query <- list( @@ -114,17 +137,21 @@ get_coverage_wcs <- function( RESPONSE_CRS = wcs_crs, ... ) + + # build and run the http request request <- build_url(url) - file <- tempfile(fileext = ".mht") - x <- GET( + mht_file <- tempfile(fileext = ".mht") + http_response <- GET( url = request, - write_disk(file) + write_disk(mht_file) ) + # multipart file extract tif part - unpack_mht(file) - file <- str_replace(file, "mht", "tif") - } + tif_file <- unpack_mht(mht_file) + } # /version 2.0.1 + + # variant: version 1.0.0 if (version == "1.0.0") { url$query <- list( SERVICE = "WCS", @@ -145,28 +172,33 @@ get_coverage_wcs <- function( RESPONSE_CRS = wcs_crs, ... ) + + # build and run the http request request <- build_url(url) - file <- tempfile(fileext = ".tif") - x <- GET( + tif_file <- tempfile(fileext = ".tif") + http_response <- GET( url = request, - write_disk(file) + write_disk(tif_file) ) } - stop_for_status(x) + # raise http errors + stop_for_status(http_response) - raster <- rast(file) + # assemble the spatial raster + raster <- rast(tif_file) template <- project(raster, output_crs) res(template) <- resolution raster <- project(raster, template) + return(raster) } -#' Unpack or extract the `tif` file part from an `mht` file +#' Unpack or extract the `tif` file part from an `mht` file. #' #' This helper function is needed on some `WCS` services from which an `mht` -#' file is downloaded rather than a `tif` file +#' file is downloaded rather than a `tif` file; returns the path to the `tif` #' #' @param path A path to the `mht` file #' @@ -174,6 +206,8 @@ get_coverage_wcs <- function( #' @importFrom assertthat assert_that #' @importFrom stringr str_detect str_replace #' +#' @return tif_path the path to the extracted geoTIFF. +#' #' @keywords internal #' #' @details Need three ways to read in the `mht` file to get the `tif` file out. @@ -195,8 +229,10 @@ unpack_mht <- function(path) { pos_end <- length(raw_vector) - (length(lines_raw[end + 1]) + 1) tif <- raw_vector[pos_start:pos_end] + tif_path <- str_replace(path, "mht", "tif") write_file( tif, - str_replace(path, "mht", "tif") + tif_path ) + return(tif_path) } diff --git a/_pkgdown.yml b/_pkgdown.yml index f724786..7bcd842 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -26,6 +26,8 @@ authors: href: "https://github.com/florisvdh" Hans Van Calster: href: "https://github.com/hansvancalster" + Falk Mielke: + href: "https://github.com/falkmielke" reference: - title: Functions for getting or displaying data from web services diff --git a/man/figures/qgis_lead.avif b/man/figures/qgis_lead.avif new file mode 100644 index 0000000..2de3cf7 Binary files /dev/null and b/man/figures/qgis_lead.avif differ diff --git a/man/get_coverage_wcs.Rd b/man/get_coverage_wcs.Rd index 95afc33..5ad289c 100644 --- a/man/get_coverage_wcs.Rd +++ b/man/get_coverage_wcs.Rd @@ -5,11 +5,11 @@ \title{Get a layer from a web coverage service within a bounding box} \usage{ get_coverage_wcs( - wcs = c("dtm", "dsm", "omz", "omw"), + wcs = c("dtm", "dsm", "omz", "omw", "dhmv"), bbox, layername, resolution, - wcs_crs = "EPSG:4258", + wcs_crs = c("EPSG:4258", "EPSG:31370"), output_crs = "EPSG:31370", bbox_crs = "EPSG:31370", version = c("1.0.0", "2.0.1"), @@ -17,7 +17,7 @@ get_coverage_wcs( ) } \arguments{ -\item{wcs}{One of \code{"dtm"}, \code{"dsm"}, \code{"omz"}, \code{"omw"}} +\item{wcs}{One of \code{"dtm"}, \code{"dsm"}, \code{"omz"}, \code{"omw"}, \code{"dhmv"}} \item{bbox}{An object of class bbox of length 4.} @@ -51,7 +51,8 @@ The following WCS services can currently be used: \item \code{"omw"}: orthophotomosaic winter images Flanders \item \code{"dtm"}: digital terrain model Flanders \item \code{"dsm"}: digital surface model Flanders -See \href{https://metadata.vlaanderen.be/metadatacenter/srv/dut/catalog.search#/search?keyword=OGC:WCS}{metadata Vlaanderen} for more information # nolint: line_length_linter. +\item \code{"dhmv"}: digital elevation model Flanders (contains dtm and dsm data) +See \href{https://metadata.vlaanderen.be/srv/eng/catalog.search#/search?any=WCS}{metadata Vlaanderen} for more information # nolint: line_length_linter. } } \examples{ diff --git a/man/unpack_mht.Rd b/man/unpack_mht.Rd index 3c594e8..c7b93f7 100644 --- a/man/unpack_mht.Rd +++ b/man/unpack_mht.Rd @@ -2,16 +2,19 @@ % Please edit documentation in R/get_coverage_wcs.R \name{unpack_mht} \alias{unpack_mht} -\title{Unpack or extract the \code{tif} file part from an \code{mht} file} +\title{Unpack or extract the \code{tif} file part from an \code{mht} file.} \usage{ unpack_mht(path) } \arguments{ \item{path}{A path to the \code{mht} file} } +\value{ +tif_path the path to the extracted geoTIFF. +} \description{ This helper function is needed on some \code{WCS} services from which an \code{mht} -file is downloaded rather than a \code{tif} file +file is downloaded rather than a \code{tif} file; returns the path to the \code{tif} } \details{ Need three ways to read in the \code{mht} file to get the \code{tif} file out. diff --git a/vignettes/spatial_dhmv_query.Rmd b/vignettes/spatial_dhmv_query.Rmd new file mode 100644 index 0000000..bc97ddf --- /dev/null +++ b/vignettes/spatial_dhmv_query.Rmd @@ -0,0 +1,502 @@ +--- +title: "Dissecting a WCS Query: DHMV Case Study." +output: rmarkdown::html_vignette +author: "Falk Mielke (falk.mielke@inbo.be)" +vignette: > + %\VignetteIndexEntry{spatial_dhmv_query} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + + + + +# DHMV Data + +Flemish government agencies host a lot of useful web services and data. Some of those are available via web services and can be found here: + +- + + +This vignette demonstrates the **query of elevation data from the Digital Elevation Model**, ``Digitaal HoogteModel Vlaanderen`` (DHMV). We will do it "the complicated way", assembling a web API query. [In the end](#orgd59526e), you will see how to use `inbospatial` to get the same result in a single function call. + +This vignette is somewhat related to, but goes beyond the use of Web Feature Services (WFS) for spatial data, [see this tutorial by Thierry Onkelinx](https://inbo.github.io/tutorials/tutorials/spatial_wfs_services). In particular, I will document some paths and workarounds to explore if errors occur in the default `ows4R` procedure. + +**TL;DR:** see the `?inbospatial:get_coverage_wcs` function to query DHMV (or other WCS) data of a given point in Flanders. An example can be found below. + + + + +# Situation/Motivation + +The purpose of the code below is to query elevation data for an arbitrary location in Flanders. It can be found here: + +- + +Unfortunately, documentation about the API is sparse, so let's hope it complies with given standards. + +We will use Web Coverage Services (WCS), as standardized by the Open Geospatial Consortium (OGC). More info on these services can be found [in the "Geocomputation with R" book, for example](https://r.geocompx.org/read-write#geographic-web-services). The go-to R package for WCS is `ows4R` by Emmanuel Blondel and contributors, [available on github](https://eblondel.github.io/ows4R). + +Unfortunately, neither the geocomputation book's example, nor the vignette in [the `ows4R` package](https://eblondel.github.io/ows4R/articles/wcs.html), nor [the INBO tutorial](https://inbo.github.io/tutorials/tutorials/spatial_wfs_services) could be transferred to the DHMV usecase. The code simply did not work. + +Therefore, I ended up going back to the nitty-gritty aspects of WCS by **manually building a query**. + + + + +# Limits of `ows4R` Package + +Before revealing [the solution](#org803f24d), I will document the error and my path to working around it: the core problem was not trivial to see, and you might run into similar issues. + +Because, in fact, the `ows4R` package seems functional. In the tutorial mentioned above, [the `ows4R` package](https://cran.r-project.org/web/packages/ows4R/index.html) is used to query data from web services. However, `ows4R` cannot handle all niche server APIs. + +First, we need to load some packages. + +```{r setup} +library("httr") +library("sf") +library("ows4R") +library("terra") +``` + +We can connect a WCS client. + +```{r} +WCS <- WCSClient$new("https://geo.api.vlaanderen.be/DHMV/wcs", "2.0.1", logger = "INFO") +``` + + +We can easily query the capabilities: + +```{r} +caps <- WCS$getCapabilities() +``` + +This would grab the same info as you can find in this xml file: + +Note that the *capabilities* xml holds a lot of valuable information: + +- First and foremost, search it for `CoverageID`; in this case, I was interested in `DHMVII_DTM_1m`: grid coverage of the LIDAR elevation data for Flanders. +- We can also find that the Coordinate Reference System (CRS) `crsSupported` is EPSG:31370, aka. [Belgian Lambert 72](https://epsg.io/31370). Always good to know. +- Finally, the capabilities xml provides a `version`: 2.0.1 at the time of writing. Yet keep that in mind for a bit. + +The "capabilities" are, per OGC standards, your go-to place for finding API documentation. + +In R, we can receive further information conveniently via `ows4R`. This works via the "coverage summary" … + +```{r} +feature <- "DHMVII_DTM_1m" +dtm_smry <- caps$findCoverageSummaryById(feature, exact = TRUE) +dtm_smry +``` + +… and the coverage description. + +```{r} +# Get description from the WCS client +dtm_des <- WCS$describeCoverage(feature) +dtm_des +``` + + +Finally, it is important to know which dimensions are available: `LON/LAT`, `x/y`, or maybe `time`? + +```{r} +dtm_dims <- dtm_smry$getDimensions() +dtm_dims +``` + + +So far, so good. + +Now, following the tutorials to get some real data: + +```{r} +x_test <- 148600 +y_test <- 208900 +range_m <- 2 # (+/- m) + +hopo_boxed <- OWSUtils$toBBOX( + xmin = x_test - range_m, + xmax = x_test + range_m, + ymin = y_test - range_m, + ymax = y_test + range_m + ) +``` + +I will write the file to disk: + +```{r, eval=FALSE} + +mht_file <- tempfile(fileext = ".mht") + +tryCatch({ + dtm_data <- dtm_smry$getCoverage( + bbox = hopo_boxed, + filename = mht_file + ) + }, error = function (err) {message(err)} +) + +raster <- rast(mht_file) + +## to solve the errors above, see below, or use: +# require("inbospatial") +# tif_file <- unpack_mht(mht_file) +# raster <- rast(tif_file) +# +# plot(raster, col = gray.colors(256)) +``` + +… or, maybe not. The function fails, somehow. + +*Status quo:* + +- I receive a file (an "mht" because *I called it* "mht"), supposedly some kind of GeoTIFF image (because I *asked* for `image/tiff`), yet wrapped into other data chunks. +- I also receive a couple of error messages and warnings: + - *"Start tag expected, '<' not found"* + - *"cannot open this file as SpatRaster:…"* + - The raw file is *"not recognized as being in a supported file format. (GDAL error 4)"* + +Maybe most useful, `ows4R` briefly displays the URL it used to attempt data query: `https://geo.api.vlaanderen.be/DHMV/wcs?service=WCS&version=2.0.1&coverageId=DHMVII_DTM_1m&subset=x(148599,148601)&subset=y(208899,208901)&format=image/tiff&request=GetCoverage` + +Let's break this URL down: + +- the desired `service` is indeed `WCS`, +- `version` is `2.0.1` (I would guess) +- `request`, `coverageID`, and `format` seem to be as expected +- there are `subset` components for the respective dimensions. + +Feel free to try yourself to paste the URL to a browser and download the file. + +Still, the downloaded file is unreadable. Subtle symptom: even if one adjusts the `subset` (i.e. `bbox` argument, above), the received tif does not change in size; the download is bbox-agnostic, so to speak. + +Opening the file with a text editor such as [vim](https://www.vim.org), we see an XML part, some generic binary stuff, some hints that this [should be a GML](https://en.wikipedia.org/wiki/Geography_Markup_Language). + +There are [sample geotiffs online](https://github.com/mommermi/geotiff_sample/blob/master/sample.tif), and they look different, or maybe not. Gml is clearly not geotiff, except maybe for the binary part (which I tried to extract, to no avail). + +The file does not download with an extension on the author's OS and browser. +On a different OS, the browser automatically extended it with ".mht". +In fact, it is [an MHT file](https://docs.fileformat.com/web/mht). +More on that [below](#mht_solution). + + + + +# The Lead: QGIS + +Well, there is advanced software to open GeoTIFF'esque files. Or geography markup. At least some software might be able to identify what we actually have here. + +First, if R failed, Python might work. + +``` +import rasterio +import rasterio.plot + +data_name = "/tmp/test_hopo.tif" +tiff = rasterio.open(data_name) +rasterio.plot.show(tiff, title = "will this work? no!") +``` + +… but it doesn't recognize the file type as raster data. + +How about **[QGIS](https://qgis.org)?** *"Spatial Without Compromise"*, they advertise. Unfortunately, that program also was unable to open the file downloaded in R. + +However, QGIS can actually connect to WCS directly. This does not help much if you need the data in R. Yet, why not. + +Look and behold: qgis can actually connect to and query elevation data for Flanders; good indication that the `geo.api` is, somehow, functional. I could zoom and move, export images. I could even export the data. I stopped exporting the data when it threatened to entirely filled my tiny system partition (the raster data is more than 100GB) - good confirmation that we want web services for this. + +And, luckily, at some careless map movement, I received an error output from QGIS. With it, QGIS gave me a working URL to query DHMV: + +- `https://geo.api.vlaanderen.be/DHMV/wcs?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&FORMAT=GeoTIFF&COVERAGE=DHMVII_DTM_1m&BBOX=162064,165735,170953,168882&CRS=EPSG:31370&RESPONSE_CRS=EPSG:31370&WIDTH=1622&HEIGHT=574` + + +```{r echo=FALSE, out.width='90%'} +knitr::include_graphics("../man/figures/qgis_lead.avif") +``` + +Let's dissect this URL. + + + + +# The Solution: Query Building + +If you start at the question mark of the URL, and split at every ampersand, you can extract the following components from the (correct) QGIS URL: + +| **correct:** | | +|---------------------- |--------------------------------- | +| SERVICE | WCS | +| VERSION | 1.0.0 | +| REQUEST | GetCoverage | +| FORMAT | GeoTIFF | +| COVERAGE | DHMVII_DTM_1m | +| BBOX | 162064,165735,170953,168882 | +| CRS | EPSG:31370 | +| RESPONSE_CRS | EPSG:31370 | +| WIDTH | 1622 | +| HEIGHT | 574 | + +Compare this to the (unsuccessful) `ows4R` attempt: + +| **non-functional:** | | +|------------------- |--------------------------------- | +| service | WCS | +| version | 2.0.1 | +| coverageId | DHMVII_DTM_1m | +| subset | x(148599,148601) | +| subset | y(208899,208901) | +| format | image/tiff | +| request | GetCoverage | + +You can selectively shuffle and adjust components, generalizing the working query and learning about the WCS API. Which feels a bit like [fuzzing](https://en.wikipedia.org/wiki/Fuzzing) the `geo.api`. + +Here is what I found: + +- `service` (`WCS`) and `request` (`DHMVII_DTM_1m`) were correct; the latter could be swapped for different datasets +- `coverage` is the correct keyword for reading `DHMVII_DTM_1m` +- `version` must be `1.0.0`; it currently does not work with 2.0.1 +- instead of subsetting, we can use a `bbox`… +- … but a `width` and `height` of the output image are mandatory; alternatively there are `resx` and `resy` (not shown) +- we must specify a `CRS` (`EPSG:31370`), optionally also a `response_crs` +- the `format` is `GeoTIFF`; slashes are weird in URLs anyway + +After some more trial and error, I could construct working URLs in R and generalize this to a function. + +```{r} + +get_elevation_wcs <- function(x, y, range_m = 1, tiff_file = NULL) { + # Please note that this function lacks a lot of + # assertions and documentation. + + # get bbox + bbox = paste(x - range_m, x + range_m, + y - range_m, y + range_m, + sep = ",") + + # the URL components + base_url = "https://geo.api.vlaanderen.be" + endpoint = "/DHMV/wcs" + + # the query parameters which worked in QGIS + elevation_query = list( + service = "WCS", + version = "1.0.0", + request = "GetCoverage", + format = "GeoTIFF", + coverage = "DHMVII_DTM_1m", + bbox = bbox, + width = 2*range_m, + height = 2*range_m, + crs = "EPSG:31370" + ) + + # get wcs data + if (is.null(tiff_file)) { + tiff_file = tempfile(fileext = ".tiff") + } + res = GET(url = modify_url(base_url, path = endpoint), + query = elevation_query, + write_disk(tiff_file)) + # note: saving this to a file is optional, + # but might prevent double download or loss of data + + # re-read raster file + data = rast(tiff_file) + + return(data) + +} +``` + +Example usage: + +```{r} +x_test <- 148600 +y_test <- 208900 +range_m <- 10 # (+/- m) +test_raster <- get_elevation_wcs(x_test, y_test, range_m = range_m) + +# we can extract a value at a point: +extract(test_raster, cbind(x_test, y_test))[[1]] + +``` + + +Voila! We can now get the elevation of a given location in Flanders via DHMV web services. + +Major grain of salt: we could not use the latest version, `v2.0.1`, of the DHMV API. +This is not a minor limitation, and should be taken seriously: old versions become obsolete, break, or stop data updates. + + + + +# Using Version `2.0.1` + +Apparently, the exact symptoms of the failing query above are somewhat depending on your browser and operating system. +Your browser might detect that the file is an `mht`, or not. + + +Here is an example of how to build a query with the current version of DHMV: + +```{r} + +url <- parse_url("https://geo.api.vlaanderen.be/DHMV/wcs") + +# the query parameters which worked in QGIS +url$query <- list( + SERVICE = "WCS", + VERSION = "2.0.1", + REQUEST = "GetCoverage", + FORMAT = "image/tiff", + COVERAGEID = "DHMVII_DTM_1m", + SUBSET = "x,http://www.opengis.net/def/crs/EPSG/0/EPSG:31370(148000,149000)", + SUBSET = "y,http://www.opengis.net/def/crs/EPSG/0/EPSG:31370(208000,209000)", + SCALEFACTOR = "1", + CRS = "EPSG:31370", + RESPONSE_CRS = "EPSG:31370" +) + +request <- build_url(url) + + +# get wcs data +mht_file <- tempfile(fileext = ".mht") + +print(request) +print(mht_file) + +response <- GET(url = request, + write_disk(mht_file)) +``` + +With the `mht` specifications [found online](https://docs.fileformat.com/web/mht), it is possible to create an unpacking function, which has become part of this very package. + +```{r} +require("readr") +require("stringr") + +unpack_mht <- function(mht_filepath) { + lines_raw <- read_lines_raw(mht_filepath) + lines_char <- suppressWarnings(read_lines(mht_filepath)) + raw_vector <- read_file_raw(mht_filepath) + + start <- which(str_detect(lines_char, "^II\\*")) + end <- length(lines_raw) - 1 + pos_start <- length(unlist(lines_raw[1:(start - 1)])) + start + pos_end <- length(raw_vector) - (length(lines_raw[end + 1]) + 1) + + tif <- raw_vector[pos_start:pos_end] + tif_path <- str_replace(mht_filepath, "mht", "tif") + write_file( + tif, + tif_path + ) + return(tif_path) +} +``` + + +This way you get the elevation data extracted: + +```{r} + +tif_file <- unpack_mht(mht_file) + +raster <- rast(tif_file) +# plot(raster, col = gray.colors(256)) + +``` + + + + + +# The "`inbospatial`" Way + +Colleague [Hans Van Calster](https://github.com/hansvancalster) had been facing the same problems with related APIs a year earlier, and he wrote a function for it in [the \`inbospatial\` package](https://github.com/inbo/inbospatial). +We now refined the `get_coverage_wcs()` procedure to solve the `mht` issue. + +Here is how to use it, for the same case as above: + +```{r} +require("sf") +require("inbospatial") + +bbox <- sf::st_bbox( + c(xmin = 148000, xmax = 149000, ymin = 208000, ymax = 209000), + crs = sf::st_crs(31370) +) +hopo_raster <- get_coverage_wcs( + wcs = "DHMV", + bbox = bbox, + layername = "DHMVII_DTM_1m", + version = "2.0.1", + wcs_crs = "EPSG:31370", + resolution = 1 +) + +plot(hopo_raster, col = gray.colors(256)) +``` + +The function gives you quite some options of APIs to query data from. + + +| **wcs** | **layer name** | **reference** | +|------------------- |--------------------------------- |---- | +| dtm | `EL.GridCoverage.DTM` | [here](https://www.vlaanderen.be/datavindplaats/catalogus/hoogte-dtm) | +| dsm | `EL.GridCoverage.DSM` | [here](https://www.vlaanderen.be/datavindplaats/catalogus/hoogte-dsm) | +| dhmv | `DHMVI_DTM_5m` | [here](https://metadata.vlaanderen.be/srv/dut/catalog.search#/metadata/9b0f82c7-57c4-463a-8918-432e41a66355) | +| dhmv | `DHMVII_DTM_1m` | [here](https://metadata.vlaanderen.be/srv/dut/catalog.search#/metadata/f52b1a13-86bc-4b64-8256-88cc0d1a8735) | +| dhmv | `DHMVII_DSM_1m` | [here](https://metadata.vlaanderen.be/srv/dut/catalog.search#/metadata/0da2e5e4-6886-426b-bb82-c0ffe6faeff6) | + +At a quick glance, we found that the `dtm` wcs seems to include canopy points, which is not according to definition; see here: A general advice is to confirm your data for a region you know. + +In addition, there are `oi-omz` and `oi-omw` for "zummer" and winter ortho-mosaic images. + +The burden of choice! + + + + +# *Post Mortem* + +The value of this vignette seems to be limited: in essence, it provides the specific query function arguments for elevation model queries with the `geo.api` of the Flemish digital services. Quite a niche use case. + +However, this is also a story of tracing errors in existing tools, and scraping other tools for hints. + +Can we actually trace the error back to a source? It was pure co-incidence that QGIS provided a working link, or are there general strategies? + +In retrospective, the initial issue was worsened by unfavorable circumstances. + +- `ows4R` retrieved a file and produced an error; another error was only thrown downstream. Neither `ows4R`, nor `terra`, `sf`, `rasterio`, or `qgis` were able to handle `mht` files. This lack of support questions the API's design choice of packing and sending an `mht`. +- Documentation for the `geo.api` is lacking (or I did not find it); in particularly I missed usage examples, vignettes, or tests. + +The colleagues at ``Digitaal Vlaanderen`` or the DHMV group certainly had enough work on their table, yet more documentation would be desirable. Same for `ows4R`, judging by their [github issue list](https://github.com/eblondel/ows4R/issues) (see e.g. [this one](https://github.com/eblondel/ows4R/issues/114)). On the other hand, [QGIS is free- and open source](https://en.wikipedia.org/wiki/QGIS), its queries seem to be well maintained and are accessible, errors are transparent, I could have read the source code - a good example piece of software. + + +Yet in general, if you run into issues like this, do not hesitate to contact support or [file a github issue](https://github.com/eblondel/ows4R/issues/114). +You should not need to fuzz an API. + + +Thank you for reading, and may all your spatial API queries be successful! + + +# Table of Contents + +- [DHMV Data](#org86b1211) +- [Situation/Motivation](#org75073fb) +- [Limits of `ows4R` Package](#orge6d66b9) +- [The Lead: QGIS](#orgfaa5a3d) +- [The Solution: Query Building](#org803f24d) +- [Using Version `2.0.1`](#mht_solution) +- [The "inbospatial" Way](#orgd59526e) +- [Post Mortem](#org575ca41) +