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

我希望它返回两件事。我的区域是否触及任何 Natura2000 区域？以及哪些区域是，换句话说，它们的名称是什么。如果我在 Qgis 中加载 WMS 层，Natura2000 区域的名称应该在 naam_n2k 下。

在我的示例中，答案应该是我的区域涉及两个 Natura2000 区域，这些 Natrua2000 区域的名称分别为 Vecht-en Beneden-Reggegebied 和 Engbertsdijksvenen。

Answer 1

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

正如@IvanSanchez 所说，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")

其次，回答您的问题：我们可以将带有几何过滤器的 CQL 查询添加到 WFS 请求中。
请注意，这在您的特定情况下是可能的，但这并不是所有 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)