Extracting data from a WMS-layer

Question

I am trying to extract data from a WMS-layer.

As an example, I would like to analyse whether or not a Natura2000 area is touched by my area and what the specifics are of the Natura2000 area.

The WMS-layer of the Natura2000 area can be found at: https://geodata.nationaalgeoregister.nl/natura2000/ows?service=WMS&request=GetLegendGraphic&format=image/png&width=20&height=20&layer=natura2000"

Let's say my area is a circle with a radius of 7500 meters around a certain x- and y-coordinate;

I was trying to get this done with the leaflet-package, but it seems more a tool to show information, instead of analysing information.

x.WGS=6.662226
y.WGS=52.53206

leaflet() %>% setView(x.WGS, y.WGS, zoom = 11) %>%   
  addTiles() %>% 
  addMarkers(lng = x.WGS, lat = y.WGS)%>%
  addWMSTiles(
    "https://geodata.nationaalgeoregister.nl/natura2000/ows?service=WMS&request=GetLegendGraphic&format=image/png&width=20&height=20&layer=natura2000",
  layers = "natura2000",
  options = WMSTileOptions(format = "image/png", transparent = TRUE),
  attribution = "") %>%
  addCircles(lng = x.WGS, lat = y.WGS, weight = 1,
  radius = 7500)

I would like it to return two things. Is any Natura2000-area touched by my area? and which areas are that, in other words, what are their names. If I load the WMS-layer in Qgis the name of an Natura2000-area should be under naam_n2k.

In my example, the answer should be that there are two Natura2000 area's touched by my area and the names of these Natrua2000 area's are Vecht- en Beneden-Reggegebied and Engbertsdijksvenen.

WMS is a raster protocol only, and the nature of it does not allow for calculating polygon intersections. I suggest you dig deeper in the datasets available at https://geodata.nationaalgeoregister.nl and see if you can locate some **vector** data. — IvanSanchez, Apr 23 '18 at 22:19

Lennert · Answer 1 · 2018-09-17T13:00:33.873

The leaflet package is a R package to access leaflet functionality, one of the most popular open-source JavaScript libraries for displaying interactive maps. More info: https://rstudio.github.io/leaflet/

As @IvanSanchez notes, a Web Map Service (WMS) is a map that you can load in your own application. See this as a georeferenced image/picture without information about the actual polygons, shapes etc.
If you want to do some analysis with polygons or retrieve information, you might need a Web Feature Service (WFS).

First of all, this is the code to load the WMS correctly.

# corrected WMS version
library(leaflet)
library(leaflet.extras) # for WMS legend

x.WGS=6.662226
y.WGS=52.53206

leaflet() %>% setView(x.WGS, y.WGS, zoom = 11) %>%
  addTiles() %>%
  addMarkers(lng = x.WGS, lat = y.WGS)%>%
  addWMSTiles(
    baseUrl = "https://geodata.nationaalgeoregister.nl/natura2000/ows?service=WMS",
    layers = "natura2000",
    options = WMSTileOptions(format = "image/png", transparent = TRUE),
    attribution = "") %>%
  addCircles(lng = x.WGS, lat = y.WGS, weight = 1,
             radius = 7500) %>%
  addWMSLegend(uri ="https://geodata.nationaalgeoregister.nl/natura2000/ows?service=WMS&request=GetLegendGraphic&format=image%2Fpng&width=20&height=20&layer=natura2000")

Secondly, to answer your question: we can add a CQL query with a geometry filter to a WFS request.
Note that this is possible in your specific case, but this is not standard for all WFS services.

# Since you want to do spatial analysis on features, we load the sf package
# more info: https://r-spatial.github.io/sf/
library(sf)

# we convert your point coordinates to the native CRS of the WFS layer:
  # make point and assign WGS84 CRS
  x <- st_point(c(x.WGS,y.WGS))
  x <- st_sfc(x) %>% st_set_crs(NA) %>% st_set_crs(4326)

  # transform to EPSG::28992 // Amersfoort RD is default CRS for the wfs layer
  x_RD <- st_transform(x,28992)

# The WFS allows spatial filtering (this is not always the case)
# so you can use a cql filter with a distance around your point:
# cql_filter=dwithin(natura2000:geom,point(241514 505696.5),7500,meters)
# combining and encoding the cql filter for a valid WFS url:
WFS_url <- paste0("http://geodata.nationaalgeoregister.nl/natura2000/wfs?",
                  "&service=wfs&version=2.0.0&request=GetFeature&",
                  "typeName=natura2000:natura2000&cql_filter=",
                  URLencode(
                    paste0("dwithin(natura2000:geom,",
                            st_as_text(x_RD),
                           ",7500, meters)"),
                    reserved = FALSE),
                  "&outputFormat=application/json"
                  )

# get WFS feature
natura2000 <- st_read(WFS_url)

# Transform to WGS84
natura2000_wgs84 <- st_transform(natura2000,4326)

# load data in leaflet
leaflet() %>%
  addTiles() %>%
  addMarkers(lng = x.WGS, lat = y.WGS) %>%
  addCircles(lng = x.WGS, lat = y.WGS, weight = 1,
             radius = 7500) %>%
  addPolygons(data = natura2000_wgs84,
              label = ~naam_n2k,
              popup = ~naam_n2k)

Extracting data from a WMS-layer

1 Answers1