WFS data importeren in R met behulp van GeoJSON

[e.j.h.polle]

Voor de Data Scientist op dit Forum:

Je kunt in R data van een WFS server inlezen met de functie functie st_read() uit het package sf.

Het advies is dan om output in (Geo)JSON formaat te vragen: outputFormat=application/json.

En vervolgens kun je deze Simple Features met de library tmap (of mapview) in een interactieve kaart bekijken.

Kijk, zo wordt GISsen met RStudio leuk

Voor meer uitleg en een voorbeeld, zie:
st_read(): WFS data importeren in R met behulp van GeoJSON

Groet,

Egge-Jan

[Jellest]

Bedankt voor het uitleggen van deze mogelijkheid.

Zelf gebruik ik de ogr2ogr functie van de gdalUtils package (in combinatie met Rgdal package) om een WFS in te lezen. Ken jij deze methode ook en zijn daar enige verschillen in qua functionaliteiten / performance?

[e.j.h.polle]

Hoi Jelle,

Dank voor je bericht. Nee, de route die jij beschrijft ken ik niet. Deze blog heb ik begin mei geschreven, en ondertussen hebben we verder uitgezocht hoe je een WFS server kunt aanroepen in R.

Zo is het handig om functie build_url() (uit de package httr) te gebruiken, zodat je de URL uit losse componenten op kan bouwen, voordat je deze aan de functie st_read() (uit de package sf) geeft.

Zie het hoofdstuk OGC Web Feature Service WFS uit onze Tutorial - Spatial Analysis in R with Open Geodata die we voor de uRos2018 conferentie geschreven hebben. Je kunt de PDF vinden in deze GitHub repository:

https://github.com/TWIAV/Spatial_Analysis_in_R_with_Open_Geodata

Er zit ook een mooi paging voorbeeld bij, zodat je in één keer meer dan 1000 records op kunt vragen.

Groet,

Egge-Jan

[Jellest]

Hoi Egge-Jan,

Nog bedankt voor je antwoord en uitleg. De tutorial is zeker zinvol en leerzaam. Ik ben bijna klaar met een BGT extract tool die via R alle GML fiies van de BGT ophaalt. Ik zal hiermee zeker de httr package gebruiken om zo zoveel mogelijk features op te halen waar nodig. ben benieuwd of het lukt.

WFS is inderdaad vooral handig voor vector data. Ik was benieuwd of je ook ervaring hebt met het aanroepen van data via een WCS in R, gezien deze thread gaat over GIS en R met PDOK data. Ik ben namelijk ook bijna klaar met een AHN extract tool die via R alle raster bladen ophaalt via atom feeds (zie AHN2 of AHN3). Daarna verkleint hij het gebied via een clip om zo maar een bepaald gebied naar keuze te downloaden. Ik vraag me af of deze route overbodig is als dit ook kan via een WCS met st_read().

En inderdaad, zo wordt GISsen met R erg leuk!

groeten,
Jelle

[Lennart]

Hoi Egge Jan,

Bedankt voor je heldere uitleg. Het werkt goed.

Heb je ook een voorbeeld met deze WFS service en de library Leaflet?
In het bovenstaande voorbeeld gebruik je tmap/mapview. Hoe werkt het met Leaflet?

Weet iemand dat toevallig? Heb je ook een voorbeeld?

Ik ben een R Shiny app met een leaflet kaartje aan het maken. Het is me gelukt om een WMS service te gebruiken in deze app. Maar het lijkt me nog toffer om ook de WFS service te gebruiken.

Groeten Lennart Jongen

[antonbakker]

Hoi Lennart, ik heb recentelijk een viewer in Leaflet gebouwd die GeoJSON ophaalt van een webservice (geen WFS, maar het principe is hetzelfde, vervang de url voor een WFS request met JSON response). De code staat hier.

[Lennart]

Hoi Anton.

Super bedankt en mooie viewer! Ik heb het in mijn R script als volgt opgelost: (ik heb de dataset Natura2000 gebieden als voorbeeld gepakt)

library(leaflet)
library(magrittr)
library(sf)

# Make WFS Url
WFS_url <- paste0("http://geodata.nationaalgeoregister.nl/natura2000/wfs?",
                  "&service=wfs&version=2.0.0&request=GetFeature",
                  "&typeName=natura2000:natura2000",
                  "&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() %>%
  addPolygons(data = natura2000_wgs84,label = ~naam_n2k,popup = ~naam_n2k)

[antonbakker]

Top! Je kan overigens ook direct de features in WGS84 opvragen bij de WFS door de volgende query parameter toe te voegen: srsName=EPSG:4326, scheelt weer aantal cpu cycles aan de client kant

[Lennart]

Bedank Anton, dat is inderdaad handig!

[Lennart]

Misschien OFF-TOPIC, maar dat leidt mij tot een tweede vraag,

Hoe kan ik de data filteren, zonder dat ik eerst de hele JSON moet downloaden?

Ik wil alleen polygonen importeren binnen een straal van 8 kilometer van een punt/locatie. Ik heb het geprobeerd met een CQL-Filter (zie script hieronder). Maar het lijkt dat de filter niet werkt en importeert de hele dataset.

# Point
x.WGS=5.696340
y.WGS=51.967951

# make point and add crs
x <- st_point(c(x.WGS,y.WGS))
x <- st_sfc(x) %>% st_set_crs(4326)

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

# Make WFS Url WITH CQL FILTER
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),
                           ",8000, meters)"),
                    reserved = FALSE)
                  ,"&srsName=EPSG:4326"
                  ,"&outputFormat=json")

# get WFS feature
natura2000_wgs84 <- st_read(WFS_url)

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

Sorry ik heb nog niet zo veel ervaring met het importeren van WFS data in R.

[e.j.h.polle]

Hoi Lennart,

Je zou eens naar onderstaande post kunnen kijken. Dit gaat niet over een straal, maar over een bounding box via WFS.

Hieronder staat de R code die ik daar als voorbeeld geef.

Tip: kijk ook eens naar de library httr (van Hadley Wickham) met de functies parse_url en build_url. Hiermee kun je al je parameters in een overzichtelijke list opgeven.

HTH,

Egge-Jan

library(httr)
library(sf)
library(tmap)

url <- parse_url("https://geodata.nationaalgeoregister.nl/nwbwegen/wfs")
url$query <- list(service = "WFS",
                  version = "2.0.0",
                  request = "GetFeature",
                  typename = "nwbwegen:wegvakken",
                  BBOX = "154500,462500,155500,463500",
                  outputFormat = "application/json")
request <- build_url(url)

nwb_wegen_amersfoort <- st_read(request)

tmap_mode("view")
tm_shape(nwb_wegen_amersfoort)+tm_lines(col="red", lwd = 4)+tm_format("NLD")

Bounding box bij WFS laag meegeven

Hallo, Ik ben een beginner mbt het inladen van gegevens via WFS in QGIS. Ik wil de wegen inladen van 1 Nederlandse stad. De wegen kan inladen door verbinding te maken met de WFS van NWB (https://geodata.nationaalgeoregister.nl/nwbwegen/wfs?request=GetCapabilities). Het is natuurlijk onnodig om alle wegen in te laden wanneer ik alleen interesse heb in de wegen van 1 stad. Is er een mogelijkheid dat ik een bounding box mee kan geven bij het opvragen van de WFS laag? En kan ik ergens documentat…

[Lennart]

Bedankt Egge Jan,

Dit werkt goed. Ik moest bij enkele datasets wel eerst kijken welk coördinaten stelsel als ‘default’ was. De BBOX coördinaten moesten bij mij in hetzelfde coördinaten stelsel staan als de default coördinaten stelsel van de betreffende dataset.

Groeten Lennart Jongen