Kaarten

[bramdhondt]

Een issue gewijd aan proper inlezen van de kaarten. Allen beschikbaar op https://www.vlaanderen.be/datavindplaats/catalogus

• [ ] Habitatrichtlijn : /habitatrichtlijndeelgebieden
• [ ] Vogelrichtlijn : /vogelrichtlijngebieden
• [ ] Habitatkaart : /biologische-waarderingskaart-en-natura-2000-habitatkaart-toestand-2020 - vervangen door toestand-2023 (onder)
• [ ] ANB-patrimoniumdatabank : /openbare-bossen-en-natuurdomeinen-beheerd-door-anb-in-naam-van-de-vlaamse-overheid
• [ ] Natuurbeheerplannen : /natuurbeheerplannen

Kan rechtstreeks ingelezen worden in R, d.m.v. de link de knop "Downloadservices:Kopieer". Bv. patrimoniumdatabank:

read_sf("https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wfs?service=WFS&version=2.0.0&request=getFeature&typeName=am:am_patdat")

[bramdhondt]

@SanderDevisscher @fleurpet : ter informatie. In de patrimoniumdatabank (ook elders?) zijn sommige geografische entiteiten van het type "MULTISURFACE" (een verzameling van punten, lijnen, polygonen...). Sommige functies kunnen daar niet mee overweg, bv. leaflet, andere weer wel, bv. ggplot.

Een workaround is om een MULTISURFACE om te zetten in een MULTPOLYGON (een verzameling van enkel polygonen). Een oplossing wordt hier gegeven.

Ik pas dat hier ter illustratie toe (uit ander project):

# patrimoniumdatabank
PatDat <- read_sf("https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wfs?service=WFS&version=2.0.0&request=getFeature&typeName=am:am_patdat")
# selectie van domein
peerd <- PatDat %>% filter(domeinnaam %in% c("Peerdsbos"))
# get crs
st_crs(peerd) # Zegt "NA"...
peerd <- st_set_crs(peerd, 31370) # ...maar is in feite Lambert (zie webpagina kaart)...
peerd <- st_transform(peerd, "+proj=longlat +datum=WGS84") # ...en moet WGS84 zijn voor leaflet.

# ggplot is geen probleem (leaflet wel)
ggplot() + geom_sf(data = peerd, fill = "lightgreen")

# van MULTISURFACE naar MULTIPOLYGON
library(gdalUtilities)
ensure_multipolygons <- function(X) {
    tmp1 <- tempfile(fileext = ".gpkg")
    tmp2 <- tempfile(fileext = ".gpkg")
    st_write(X, tmp1)
    ogr2ogr(tmp1, tmp2, f = "GPKG", nlt = "MULTIPOLYGON")
    Y <- st_read(tmp2)
    st_sf(st_drop_geometry(X), geom = st_geometry(Y))
}
peerd <- ensure_multipolygons(peerd)

# leaflet lukt nu wel
leaflet() %>% 
  setView(lat = 51.273568, lng = 4.485597, zoom = 10) %>%
  addTiles(group = "OSM") %>%
  addPolygons(data = peerd)

[fleurpet]

@bramdhondt: Ik zag dat er sinds vorige week een vernieuwde versie van de BWK en Natura 2000 kaart online staat: Biologische Waarderingskaart en Natura 2000 Habitatkaart - Toestand 2023. Gebruik ik deze i.p.v. die van 2020?

[bramdhondt]

@bramdhondt: Ik zag dat er sinds vorige week een vernieuwde versie van de BWK en Natura 2000 kaart online staat: Biologische Waarderingskaart en Natura 2000 Habitatkaart - Toestand 2023. Gebruik ik deze i.p.v. die van 2020?

@fleurpet - Ja, lijkt me goed!

[fleurpet]

Ik heb alle WFS kaarten kunnen inladen.

Alleen lijkt de output voor de BWK en Natura 2000 habitatkaart niet te kloppen.

Ik heb de WFS ingeladen volgens volgend script:

wfs_bwk <- "https://geo.api.vlaanderen.be/BWK/wfs"

url <- parse_url(wfs_bwk)
url$query <- list(service = "wfs",
                  version = "2.0.0", # facultative
                  request = "GetFeature",
                  typename = "BWK:Bwkhab")
request <- build_url(url)
request

bwkhab <- read_sf(request) %>%
  st_transform("+proj=longlat +datum=WGS84") %>%  
  ensure_multipolygons()   

vb. voor zone rond Gent:

Dit is mijn resultaat in leaflet (habitatkaart in paars):

Dit is de BWK en Natura 2000 Habitatkaart van Geopunt:

@bramdhondt @SanderDevisscher - Weet iemand van jullie waar het fout gaat?

[bramdhondt]

@fleurpet -- Inderdaad, ik bekom hetzelfde beeld, en wel via onafhankelijke weg. FYI, ik deed het wat anders:

# Link van www.vlaanderen.be-pagina ("Downloadservices" -> "Kopieer")
# MAAR: "&maxFeatures=1" uit de URL gehaald
x <- read_sf("https://geo.api.vlaanderen.be/BWK/wfs?service=WFS&version=1.1.0&request=GetFeature&typeName=BWK:Bwkhab&outputFormat=application/json&srsName=epsg:3857")
y <- st_set_crs(x, 31370)
y <- st_transform(x, "+proj=longlat +datum=WGS84")
# geen MULTISURFACE aanwezig, dus "ensure_multipolygons(y)" niet strikt nodig
leaflet() %>% 
  setView(lat = 51.273568, lng = 4.485597, zoom = 6) %>%
  addTiles(group = "OSM") %>%
  addPolygons(data = y)

Dit is volgens mij het beeld dat we moeten bekomen (geopend via Geopunt). (Habitatcodes zijn maar voor een kleiner aantal percelen gedefinieerd dan BWK-codes.)

Wat mij heel erg opvalt, is dat de kaartlaag exact 10.000 polygonen bevat. What are the odds? @fleurpet , is dit bij jou ook het geval? @SanderDevisscher , enig idee of er hier iets van een plafonnering gebeurt? Anders zou ik dat eens aan Carine Wils voorleggen.

[fleurpet]

@bramdhondt - Het viel me inderdaad ook op dat het exact 10.000 polygonen zijn. Ik heb de kaart eens geopend in ArcMap en zie dat de oorspronkelijke shapefile 70.826 polygonen bevat. Het lijkt er dus inderdaad op dat er een cut-off staat op 10.000.

[bramdhondt]

@fleurpet - Goed opgehelderd! Ga jij op zoek naar een oplossing (bv. online, teamgenoten, Carine Wils, Floris Vanderhaeghe ...)?

[SanderDevisscher]

Er staat blijkbaar een limiet op het aantal features dat je via R vanuit een wfs kan inlezen (en ja dat is 10000). Ik zie niet direct een oplossing muz de shape aanmaken via arcgis en die ergens parkeren en inlezen.

[fleurpet]

Het is gelukt om het in te lezen als shapefile i.p.v. WFS. Nu klopt de output:

Aangezien deze shapefile 70.826 features bevat is deze shapefile wel best zwaar en lijkt me niet heel gebruiksvriendelijk. @SanderDevisscher Is er een manier om deze makkelijker te maken in gebruik? Eventueel verschillende polygonen samen nemen?

[SanderDevisscher]

@fleurpet ik ken de structuur van de shape niet zo goed maar als er een kolom is waarop we de polygonen kunnen groeperen, kunnen we deze idd dissolven (zie https://www.jla-data.net/2018/02/14/2018-02-14-dissolving-polygons-in-sf-environment/ voor een voorbeeld werkwijze)