Als ik op een perceel klik haal ik via een WFS request de bijbehorende data op. Nu wil ik graag nog een los plaatje ophalen van het zojuist aangeklikte perceel om deze als image te plaatsen. Ik heb gezien dat dit via een WMS request kan op deze manier:
Alleen hoe vind ik nu na een klik op een perceel (in dit geval perceel “266”) de bijbehorende BBOX coordinaten om deze request te voltooien?
Kun je daar de data die je uit de WFS haalt niet voor gebruiken? Kwestie van de min/max x en y van de geometrie vinden.
Beste Simeon,
Hoe herleidt je dan de min/max x,y waarden, uit de opgehaalde coördinaten array?
ik ben al aan het zoeken geweest naar een VBA-script (voor gebruik in excel), want daar gebruik ik een XMLFILTERWebservce( GetFeature-url ://geom) op de feature, omdat ik geen rechtstreekse functie gevonden heb.
Hoi Robert,
Het is al heel lang geleden dat ik iets in VBA heb gedaan, dus ik beschijf het even:
Je hebt al een array. Neem de 1e x en y, en zet die in een maxX, maxY, minX en minY variabele. Loop door je array, en als de betreffende waarden respectievelijk groter of kleiner zijn, zet je ze in je min- en max-variabelen. Als je door je hele array heen bent gelopen, bevatten je min- en max-variabelen de grootste en kleinste waardes, en dat is dus je MBR (minimum bounding rectangle).
Overigens ben ik wel verbaasd dat je dit in Excel wil doen, QGis (of een ander gis- of zelfs cad-pakket) zou dit ook, en m.i. nog veel makkelijker kunnen doen. Afhankelijk van wat je wil bereiken, vermoed ik dat QGis het out of the box kan (eventueel met behulp van de PDOK-plugin). [Maar als je dat zou willen proberen moeten we maar even een apart draadje starten.
De GML output van de WFS geeft een boundingbox terug in de response, dus dat scheelt een bbox afleiden van een feature.
Zie hieronder een Python voorbeeld script voor het genereren van een GetMap url op basis van een WFS GetFeature response.
Input:
https://geodata.nationaalgeoregister.nl/kadastralekaart/wfs/v4_0?request=GetFeature&service=WFS&typenames=kadastralekaartv4:perceel&version=2.0.0&filter=%3CFilter%3E%3CPropertyIsEqualTo%3E%3CPropertyName%3EidentificatieLokaalID%3C/PropertyName%3E%3CLiteral%3E13040217570000%3C/Literal%3E%3C/PropertyIsEqualTo%3E%3C/Filter%3E
Geeft output:
https://geodata.nationaalgeoregister.nl/kadastralekaart/wms/v4_0?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=false&LAYERS=Perceelvlak,label&CRS=EPSG%3A28992&WIDTH=400&HEIGHT=495.84326643824795&BBOX=140986.204,470626.374,141009.275,470654.973
#!/usr/bin/env python3
import requests
from lxml import etree
lokaal_id="13040217570000"
getfeature_url=f"https://geodata.nationaalgeoregister.nl/kadastralekaart/wfs/v4_0?request=GetFeature&service=WFS&typenames=kadastralekaartv4:perceel&version=2.0.0&filter=%3CFilter%3E%3CPropertyIsEqualTo%3E%3CPropertyName%3EidentificatieLokaalID%3C/PropertyName%3E%3CLiteral%3E{lokaal_id}%3C/Literal%3E%3C/PropertyIsEqualTo%3E%3C/Filter%3E"
print(f"GetFeature URL: {getfeature_url}")
response = requests.get(getfeature_url)
nsmap = {'gml': 'http://www.opengis.net/gml/3.2',"wfs":"http://www.opengis.net/wfs/2.0"}
root = etree.XML(response.content)
lc = next(iter(root.xpath('//wfs:boundedBy/gml:Envelope/gml:lowerCorner/text()', namespaces=nsmap)),None) # xpath returns list, get first item or None if empty
uc = next(iter(root.xpath('//wfs:boundedBy/gml:Envelope/gml:upperCorner/text()', namespaces=nsmap)),None)
bbox_string = [*lc.split(" "), *uc.split(" ")] # unpack lists and combine in new list with * operator
bbox=[float(x) for x in bbox_string]
bbox_pad = [bbox[0]-5,bbox[1]-5,bbox[2]+5,bbox[3]+5] # add 5 meter padding to bounds
d_x = bbox_pad[2] - bbox_pad[0]
d_y = bbox_pad[3] - bbox_pad[1]
width=400 # restrict width getmap request to 400px, and calculate height based on that
height=(width/d_x)*d_y
bbox_pad_string = ",".join([str(x) for x in bbox_pad])
getmap_url = f"https://geodata.nationaalgeoregister.nl/kadastralekaart/wms/v4_0?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=false&LAYERS=Perceelvlak,label&CRS=EPSG%3A28992&WIDTH={width}&HEIGHT={height}&BBOX={bbox_pad_string}"
print(f"GetMap URL: {getmap_url}")
Ha! Da’s waar ook! Ik werk eigenlijk alleen nog maar met geojson, omdat ik dat makkelijker vind. En daar zit geen bounding box in, ik was straal vergeten dat dat standaard in de GML-output zit.
Dat scheelt inderdaad een hoop werk…
Zit inderdaad wel weer meer werk in het parsen van de XML t.o.v. GeoJSON. Dus dat is een afweging om te maken.
Er staat een ticket op de sprint van tweede helft van mei om bounding box toe te voegen aan de WFS JSON output. Kan zijn dat er nog geschoven wordt in de planning. We zullen hier een berichtje posten als de bounding box is toegevoegd aan de WFS JSON output.
Zojuist uitgerold! De services op service.pdok.nl bevatten vanaf heden in de geojson ook een bbox.