从WMS层中提取数据

Question

我正在尝试从WMS层中提取数据。

例如，我想分析一下Natura2000区域是否被我的区域所触及，以及Natura2000区域的具体内容。

Natura2000区域的WMS层可以在：https://geodata.nationaalgeoregister.nl/natura2000/ows?service=WMS&request=GetLegendGraphic&format=image/png&width=20&height=20&layer=natura2000中找到。

假设我的面积是一个半径为7500米的圆，围绕着某个x和y坐标；

我试着用传单来完成这件事，但它似乎更像是一种显示信息的工具，而不是分析信息。

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)

我想要它还两件东西。有任何自然保护区被我的区域接触过吗？换句话说，它们的名字是什么？如果我在Qgis中加载WMS层，那么Natura2000 2000区域的名称应该在naam_n2k下。

在我的例子中，答案应该是有两个Natura2000区域被我的区域所触及，而这些Natrua2000区域的名称是Vecht Beneden-Reggegebied和

Answer 1

传单包是一个访问传单功能的R包，它是用于显示交互式地图的最流行的开源JavaScript库之一。更多信息：https://rstudio.github.io/leaflet/

正如@IvanSanchez所指出的，Web地图服务(Web，WMS)是您可以在自己的应用程序中加载的映射。将此视为地理参考图像/图片，而不提供实际多边形、形状等信息。

如果要使用多边形进行分析或检索信息，则可能需要Web功能服务(WFS)。

首先，这是正确加载WMS的代码。

# 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")

其次，要回答您的问题:我们可以向WFS请求添加带有几何过滤器的CQL查询。

请注意，这在您的特定情况下是可能的，但这并不是所有WFS服务的标准。

# 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)