If you have worked with digital elevation models (dem), you probably downloaded big files and had to cut out your region of interest. Or you stitched together multiple files to get the desired coverage. With Web Coverage Service (WCS) you can download raster files in exactly the extent and resolution you want to work with.

We will show how to use the wcs service of the German state Nordrhein-Westfalen (NRW). If you can read German you can read more in the manual about wcs.

We will be working with a digtial elevation model (dem or dgm in German). We will also be using digital ortho photos (dop).

get a bounding box

In order to download data we need a bounding box. We can use the the osmdata package to search for coordinates of a landmark (or region) with the getbb() function. We will be searching for the “Drachenfels”, a famous mountaintop near the river Rhine.

library(osmdata)
library(units)
library(sf)

WGS_crs <- st_crs(4326)    # coordinate system of the bbox function
UTM_crs <- st_crs(25832)   # coordinate system of the dem query

# download bounding box as matrix
Drachenfels <- getbb("Drachenfels", featuretype="landmark", format_out = "matrix", limit=1)

Drachenfels_1km_UTM <- 
  Drachenfels[,"min"] %>%               # use min point as "center point"
  st_point() %>%                        # convert to point object
  st_sfc() %>%                          # encapsulate in sfc (simple feature geometry) object
  st_set_crs(WGS_crs) %>%               # give it a crs (coordinate reference system)
  st_transform(UTM_crs) %>%             # convert to target crs
  st_buffer(dist=set_units(1, km)) %>%  # add a buffer of 1km around it
  st_bbox() %>%                         # calculate bounding box of buffered zone
  round()       # round coordinates (important, since the data query doesn't work with floats)

What we did here is get the coordinates with getbb() as a matrix of two points. Then we convert the min point to an st_point object and further convert it to an sfc object from the sf package. We then have to give it a crs (since the getbb function doesn’t add this) with st_set_crs. The query we will use later requires an UTM coordinate system. So we transform to that crs with st_transform. We then add a buffer around this central coordinate with st_buffer and calculate the bounding box of that with st_bbox.

build the query url

Now that we have a bounding box we can insert its coordinates into the data query url using the str_glue() function of the stringr package. str_glue() replaces everything inside curly braces.

library(stringr)

url <- str_glue("https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?VERSION=2.0.1&SERVICE=wcs&REQUEST=GetCoverage&COVERAGEID=nw_dgm&FORMAT=image/tiff&SUBSET=x({Drachenfels_1km_UTM$xmin},{Drachenfels_1km_UTM$xmax})&SUBSET=y({Drachenfels_1km_UTM$ymin},{Drachenfels_1km_UTM$ymax})&SCALEFACTOR=1&SUBSETTINGCRS=EPSG:25832")

The resulting url is the following: https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?VERSION=2.0.1&SERVICE=wcs&REQUEST=GetCoverage&COVERAGEID=nw_dgm&FORMAT=image/tiff&SUBSET=x(372511,374511)&SUBSET=y(5613116,5615116)&SCALEFACTOR=1&SUBSETTINGCRS=EPSG:25832. You can also use it in the browser to download the file.

The wcs url consist of the following parts:

The base url: https://www.wcs.nrw.de/geobasis/wcs_nw_dgm
a ? followed by all the options (devided by &)
VERSION=2.0.1 defines the version of the protocol
SERVICE=wcs defines that we want to use the wcs service
REQUEST=GetCoverage is the request to get data. Other requests are for example GetCapabilities or DescribeCoverage
COVERAGEID=nw_dgm defines the dataset we want to download
FORMAT=image/tiff defines the output format
SUBSET=x(372511,374511) and SUBSET=y(5613116,5615116) define the bounding box from which we want to get data
SUBSETTINGCRS=EPSG:25832 defines the coordinate system of the subset bounding box coordinates. Some services support multiple coordinate systems. Others work only with one specific.
SCALEFACTOR=1 can be used to scale down data before the download. Use SCALEFACTOR=0.1 for 10x aggregation for example.

Further options not used here are:

OUTPUTCRS=EPSG:25832 defines the output crs. For supported crs see the GetCapabilities request below.
INTERPOLATION=nearest defines the interpolation method for scaled downloads. This can also be bilinear or average or other values documentend in GetCapabilities.
RANGESUBSET=band1,band2 defines the bands to download (relevant for the orthophotos).

More options are described in the official manual on page 8.

download and read raster data

We can now download the data and read it with the raster package

library(raster)

dem <- raster(url)

names(dem) <- "elevation"

With names(dem) we give the layer a sensible name.

plot(dem)

plot of chunk plot dem

GetCapabilities

To get more information about the data service we can use the GetCapabilities request. This returns a XML file that we can read like this:

library(xml2)
capabilities <- read_xml("https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?VERSION=2.0.1&SERVICE=WCS&REQUEST=GetCapabilities")

We can use the htmltidy package to browse through the xml file:

#remotes::install_git("https://git.rud.is/hrbrmstr/htmltidy.git")
library(htmltidy)
xml_view(capabilities)

<?xml version="1.0" encoding="UTF-8"?>
<wcs:Capabilities xmlns:wcs="http://www.opengis.net/wcs/2.0" xmlns:ows="http://www.opengis.net/ows/2.0" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:gml="http://www.opengis.net/gml/3.2" xmlns:gmlcov="http://www.opengis.net/gmlcov/1.0" xmlns:swe="http://www.opengis.net/swe/2.0" xmlns:inspire_common="http://inspire.ec.europa.eu/schemas/common/1.0" xmlns:inspire_dls="http://inspire.ec.europa.eu/schemas/inspire_dls/1.0" xmlns:crs="http://www.opengis.net/wcs/crs/1.0" xmlns:int="http://www.opengis.net/wcs/interpolation/1.0" xsi:schemaLocation="http://www.opengis.net/wcs/2.0 http://schemas.opengis.net/wcs/2.0/wcsAll.xsd http://inspire.ec.europa.eu/schemas/inspire_dls/1.0 http://inspire.ec.europa.eu/schemas/inspire_dls/1.0/inspire_dls.xsd" version="2.0.1">
  <ows:ServiceIdentification>
    <ows:Title>WCS NW DGM</ows:Title>
    <ows:Abstract>Höhenmodell des Landes NRW.</ows:Abstract>
    <ows:Keywords>
      <ows:Keyword>NW</ows:Keyword>
      <ows:Keyword>NRW</ows:Keyword>
      <ows:Keyword>Nordrhein-Westfalen</ows:Keyword>
      <ows:Keyword>Bezirksregierung Köln</ows:Keyword>
      <ows:Keyword>Abteilung 7</ows:Keyword>
      <ows:Keyword>Geobasis NRW</ows:Keyword>
      <ows:Keyword>Geobasisdaten</ows:Keyword>
      <ows:Keyword>Landesvermessung</ows:Keyword>
      <ows:Keyword>AdV</ows:Keyword>
      <ows:Keyword>Arbeitsgemeinschaft der Vermessungsverwaltungen der Länder</ows:Keyword>
      <ows:Keyword>AdV-OWS-Basisprofil</ows:Keyword>
      <ows:Keyword> AdV-WCS-Profil 2.0.1</ows:Keyword>
      <ows:Keyword>WCS</ows:Keyword>
      <ows:Keyword>WCS_NW_DGM</ows:Keyword>
      <ows:Keyword>DGM</ows:Keyword>
      <ows:Keyword>Digitales Geländemodell</ows:Keyword>
      <ows:Keyword>Geländemodell</ows:Keyword>
      <ows:Keyword>INSPIRE</ows:Keyword>
    </ows:Keywords>
    <ows:ServiceType codeSpace="OGC">OGC WCS</ows:ServiceType>
    <ows:ServiceTypeVersion>2.0.1</ows:ServiceTypeVersion>
    <ows:ServiceTypeVersion>1.1.1</ows:ServiceTypeVersion>
    <ows:ServiceTypeVersion>1.0.0</ows:ServiceTypeVersion>
    <ows:Profile>http://www.opengis.net/spec/WCS/2.0/conf/core</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/WCS_protocol-binding_get-kvp/1.0/conf/get-kvp</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/WCS_protocol-binding_post-xml/1.0/conf/post-xml</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/GMLCOV/1.0/conf/gml-coverage</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/GMLCOV/1.0/conf/multipart</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/GMLCOV/1.0/conf/special-format</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/GMLCOV_geotiff-coverages/1.0/conf/geotiff-coverage</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/WCS_service-extension_crs/1.0/conf/crs</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/WCS_service-extension_scaling/1.0/conf/scaling</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/WCS_service-extension_range-subsetting/1.0/conf/record-subsetting</ows:Profile>
    <ows:Profile>http://www.opengis.net/spec/WCS_service-extension_interpolation/1.0/conf/interpolation</ows:Profile>
    <ows:Fees>Nutzungsbedingungen: Die Geobasisdaten des amtlichen Vermessungswesens werden als öffentliche Aufgabe gem. VermKatG NRW und gebührenfrei nach Open Data-Prinzipien über online-Verfahren bereitgestellt. Nutzungsbedingungen: siehe https://www.bezreg-koeln.nrw.de/system/files/media/document/file/lizenzbedingungen_geobasis_nrw.pdf</ows:Fees>
    <ows:AccessConstraints>NONE</ows:AccessConstraints>
  </ows:ServiceIdentification>
  <ows:ServiceProvider>
    <ows:ProviderName>Geobasis NRW</ows:ProviderName>
    <ows:ProviderSite xlink:type="simple" xlink:href="http://www.geobasis.nrw.de"/>
    <ows:ServiceContact>
      <ows:IndividualName></ows:IndividualName>
      <ows:PositionName></ows:PositionName>
      <ows:ContactInfo>
        <ows:Phone>
          <ows:Voice>+49(0)221-147-4994</ows:Voice>
          <ows:Facsimile>+49(0)221-147-4874</ows:Facsimile>
        </ows:Phone>
        <ows:Address>
          <ows:DeliveryPoint>Muffendorfer Str. 19-21</ows:DeliveryPoint>
          <ows:City>Bonn</ows:City>
          <ows:AdministrativeArea>Nordrhein-Westfalen</ows:AdministrativeArea>
          <ows:PostalCode>53177</ows:PostalCode>
          <ows:Country>Deutschland</ows:Country>
          <ows:ElectronicMailAddress>geobasis@bezreg-koeln.nrw.de</ows:ElectronicMailAddress>
        </ows:Address>
        <ows:OnlineResource xlink:type="simple" xlink:href="http://www.geobasis.nrw.de"/>
        <ows:HoursOfService>Montag bis Donnerstag 8:30 - 15:00</ows:HoursOfService>
        <ows:ContactInstructions></ows:ContactInstructions>
      </ows:ContactInfo>
      <ows:Role>Ansprechpartner</ows:Role>
    </ows:ServiceContact>
  </ows:ServiceProvider>
  <ows:OperationsMetadata>
    <ows:Operation name="GetCapabilities">
      <ows:DCP>
        <ows:HTTP>
          <ows:Get xlink:type="simple" xlink:href="https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?"/>
          <ows:Post xlink:type="simple" xlink:href="https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?">
            <ows:Constraint name="PostEncoding">
              <ows:AllowedValues>
                <ows:Value>XML</ows:Value>
              </ows:AllowedValues>
            </ows:Constraint>
          </ows:Post>
        </ows:HTTP>
      </ows:DCP>
    </ows:Operation>
    <ows:Operation name="DescribeCoverage">
      <ows:DCP>
        <ows:HTTP>
          <ows:Get xlink:type="simple" xlink:href="https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?"/>
          <ows:Post xlink:type="simple" xlink:href="https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?">
            <ows:Constraint name="PostEncoding">
              <ows:AllowedValues>
                <ows:Value>XML</ows:Value>
              </ows:AllowedValues>
            </ows:Constraint>
          </ows:Post>
        </ows:HTTP>
      </ows:DCP>
    </ows:Operation>
    <ows:Operation name="GetCoverage">
      <ows:DCP>
        <ows:HTTP>
          <ows:Get xlink:type="simple" xlink:href="https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?"/>
          <ows:Post xlink:type="simple" xlink:href="https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?">
            <ows:Constraint name="PostEncoding">
              <ows:AllowedValues>
                <ows:Value>XML</ows:Value>
              </ows:AllowedValues>
            </ows:Constraint>
          </ows:Post>
        </ows:HTTP>
      </ows:DCP>
    </ows:Operation>
    <ows:ExtendedCapabilities>
      <inspire_dls:ExtendedCapabilities>
        <inspire_common:MetadataUrl xsi:type="inspire_common:resourceLocatorType">
          <inspire_common:URL>https://apps.geoportal.nrw.de/soapServices/CSWStartup?Service=CSW&amp;Request=GetRecordById&amp;Version=2.0.2&amp;outputSchema=http://www.isotc211.org/2005/gmd&amp;elementSetName=full&amp;id=9d9ffbae-3f27-437a-92cd-5550af618690</inspire_common:URL>
          <inspire_common:MediaType>application/xml</inspire_common:MediaType>
        </inspire_common:MetadataUrl>
        <inspire_common:SupportedLanguages>
          <inspire_common:DefaultLanguage>
            <inspire_common:Language>ger</inspire_common:Language>
          </inspire_common:DefaultLanguage>
        </inspire_common:SupportedLanguages>
        <inspire_common:ResponseLanguage>
          <inspire_common:Language>ger</inspire_common:Language>
        </inspire_common:ResponseLanguage>
        <inspire_dls:SpatialDataSetIdentifier>
          <inspire_common:Code>https://registry.gdi-de.org/id/de.nw/DGM1</inspire_common:Code>
        </inspire_dls:SpatialDataSetIdentifier>
      </inspire_dls:ExtendedCapabilities>
    </ows:ExtendedCapabilities>
  </ows:OperationsMetadata>
  <wcs:ServiceMetadata>
    <wcs:formatSupported>image/tiff</wcs:formatSupported>
    <wcs:formatSupported>image/png</wcs:formatSupported>
    <wcs:formatSupported>image/jpeg</wcs:formatSupported>
    <wcs:formatSupported>image/png; mode=8bit</wcs:formatSupported>
    <wcs:formatSupported>image/vnd.jpeg-png</wcs:formatSupported>
    <wcs:formatSupported>image/vnd.jpeg-png8</wcs:formatSupported>
    <wcs:Extension>
      <int:InterpolationMetadata>
        <int:InterpolationSupported>NEAREST</int:InterpolationSupported>
        <int:InterpolationSupported>AVERAGE</int:InterpolationSupported>
        <int:InterpolationSupported>BILINEAR</int:InterpolationSupported>
      </int:InterpolationMetadata>
      <crs:CrsMetadata>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/3034</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/3035</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/3043</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/3044</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/3045</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/3857</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/4258</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/4326</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/4647</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/5649</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/5650</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/5651</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/5652</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/5653</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/28992</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/25831</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/25832</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/25833</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/31466</crs:crsSupported>
        <crs:crsSupported>http://www.opengis.net/def/crs/EPSG/0/31467</crs:crsSupported>
      </crs:CrsMetadata>
    </wcs:Extension>
  </wcs:ServiceMetadata>
  <wcs:Contents>
    <wcs:CoverageSummary>
      <wcs:CoverageId>nw_dgm</wcs:CoverageId>
      <wcs:CoverageSubtype>RectifiedGridCoverage</wcs:CoverageSubtype>
      <ows:Metadata xlink:type="simple" xlink:href="https://apps.geoportal.nrw.de/soapServices/CSWStartup?Service=CSW&amp;Request=GetRecordById&amp;Version=2.0.2&amp;outputSchema=http://www.isotc211.org/2005/gmd&amp;elementSetName=full&amp;id=0c6796e5-9eca-4ae6-8b32-1fcc5ae5c481" about="TC211" xlink:role="application/xml"/>
    </wcs:CoverageSummary>
  </wcs:Contents>
</wcs:Capabilities>

We can for example extract a list of all the supported crs from the xml file:

capabilities %>%
  xml_find_all(".//crs:crsSupported") %>%
  xml_text %>% cat(sep = "\n")

## http://www.opengis.net/def/crs/EPSG/0/3034
## http://www.opengis.net/def/crs/EPSG/0/3035
## http://www.opengis.net/def/crs/EPSG/0/3043
## http://www.opengis.net/def/crs/EPSG/0/3044
## http://www.opengis.net/def/crs/EPSG/0/3045
## http://www.opengis.net/def/crs/EPSG/0/3857
## http://www.opengis.net/def/crs/EPSG/0/4258
## http://www.opengis.net/def/crs/EPSG/0/4326
## http://www.opengis.net/def/crs/EPSG/0/4647
## http://www.opengis.net/def/crs/EPSG/0/5649
## http://www.opengis.net/def/crs/EPSG/0/5650
## http://www.opengis.net/def/crs/EPSG/0/5651
## http://www.opengis.net/def/crs/EPSG/0/5652
## http://www.opengis.net/def/crs/EPSG/0/5653
## http://www.opengis.net/def/crs/EPSG/0/28992
## http://www.opengis.net/def/crs/EPSG/0/25831
## http://www.opengis.net/def/crs/EPSG/0/25832
## http://www.opengis.net/def/crs/EPSG/0/25833
## http://www.opengis.net/def/crs/EPSG/0/31466
## http://www.opengis.net/def/crs/EPSG/0/31467

Or the supported interpolation methods:

capabilities %>%
  xml_find_all(".//int:InterpolationSupported") %>%
  xml_text %>% cat(sep = "\n")

## NEAREST
## AVERAGE
## BILINEAR

Or the output formats supported:

capabilities %>%
  xml_find_all(".//wcs:formatSupported") %>%
  xml_text %>%cat(sep = "\n")

## image/tiff
## image/png
## image/jpeg
## image/png; mode=8bit
## image/vnd.jpeg-png
## image/vnd.jpeg-png8

But note: The documentation states that only tiff files are geo referenced.

Finally we can get the coverage ID that we need to set in the query:

capabilities %>%
  xml_find_all(".//wcs:CoverageId") %>%
  xml_text %>% cat(sep = "\n")

## nw_dgm

download digital orthophotos (dop)

Since Orthophotos are multichannel files (for red, green, blue and infrared) we have to use the stack() function instead of raster() to read them. The query is similar to the dem. We add RANGESUBSET=1,2,3 to the query to omit the infrared band (band 4). This saves some download size.

dop <- stack(str_glue("https://www.wcs.nrw.de/geobasis/wcs_nw_dop?VERSION=2.0.1&SERVICE=wcs&REQUEST=GetCoverage&COVERAGEID=nw_dop&FORMAT=image/tiff&SUBSET=x({Drachenfels_1km_UTM$xmin},{Drachenfels_1km_UTM$xmax})&SUBSET=y({Drachenfels_1km_UTM$ymin},{Drachenfels_1km_UTM$ymax})&SCALEFACTOR=0.1&RANGESUBSET=1,2,3"))

names(dop) <- c("r","g","b")

Since the dop has 10 times the resolution of the dem we set SCALEFACTOR=0.1. This way we get a raster with the same pixel size as the dem.

Similar to the dem we add names to the layers (red, green and blue).

If we look into GetCapabilities of the dop service we can see that this one only supports one crs:

dop_capabilities <- read_xml("https://www.wcs.nrw.de/geobasis/wcs_nw_dop?VERSION=2.0.1&SERVICE=WCS&REQUEST=GetCapabilities") 

dop_capabilities %>%
  xml_find_all(".//crs:crsSupported") %>%
  xml_text %>% cat(sep = "\n")

## http://www.opengis.net/def/crs/EPSG/0/25832

plot 3d graphic with orthophoto as texture

The rayshader package is a wonderful package to create a 3d plot from a dem. You have to install it from github via devtools::install_github("tylermorganwall/rayshader") since it is not yet on CRAN.

library(rayshader)

In order to plot the dem with the orthophotos as texture we have to do some conversion as described in this tutorial:

dop_r_matrix = raster_to_matrix(dop$r)
dop_g_matrix = raster_to_matrix(dop$g)
dop_b_matrix = raster_to_matrix(dop$b)
dem_matrix = raster_to_matrix(dem)

siebeng_rgb_array = array(0,dim=c(nrow(dop_r_matrix),ncol(dop_r_matrix),3))

siebeng_rgb_array[,,1] = dop_r_matrix/255 #Red layer
siebeng_rgb_array[,,2] = dop_g_matrix/255 #Blue layer
siebeng_rgb_array[,,3] = dop_b_matrix/255 #Green layer

siebeng_rgb_array = aperm(siebeng_rgb_array, c(2,1,3)) # rotate image to correct orientation

Now we can create an interactive rotatable 3d-plot with the plot_3d function:

plot_3d(siebeng_rgb_array, dem_matrix, zscale = 1, theta = -45, phi=15, fov = 50, triangulate=T, max_error=1, windowsize=c(1400, 800), zoom=0.5)

render_snapshot()

plot of chunk snapshot 3d scene

To better the performance you can set triangulate=T and max_error to a higher number. This way the graphic will consist of fewer triangles and hence will be faster to render. But be aware that you loose some details with really high values of max_error:

plot_3d(siebeng_rgb_array, dem_matrix, zscale = 1, theta = -45, phi=15, fov = 50, zoom=.5,windowsize=c(1400, 800), triangulate=T, max_error=100)

plot of chunk snapshot 3d scene lower resolution

We can also render a movie:

plot_3d(siebeng_rgb_array, dem_matrix, zscale = 1, theta = 45, phi=30, fov = 50, triangulate=T, max_error=1, windowsize=c(1400, 800), zoom=0.5)
render_movie(filename = "../figure/source/2023-07-10-wcs-data/Drachenfels_orbit.gif", phi=20, zoom = 0.5, width=1400, height=1000)