Ik ben op zoek naar een manier om in python de hoogte van een gebouw te bepalen op een specifieke positie (het liefst dmv coordinaten). Via PDOK 3D hoogtestatistieken lukt het me om de hoogte van gebouwen te laden alleen staan hier geen locaties bij. Weet iemand hoe ik hier de locaties van kan vinden of is er misschien een andere manier waarop ik aan hoogte data per locatie in Nederland kan komen?
Je kan met de WCS service het AHN voor elke locatie in Nederland aanroepen en daar zelf eventueel de hoogte van panden uit afleiden.
Wat achtergrond: Starten met AHN WCS service - AHN - Geoforum
Ha @Tjom96, wat betreft de 3D hoogstestatitieken zou je eenvoudig de GeoPackage moeten kunnen bevragen met een SQL query, ware het niet dat de GeoPackage die ter beschikking is geen geldig SQLite bestand is. Omdat de identifier die gebruikt wordt voor de tabelnaam niet geldig is, deze start namelijk met een numeriek teken: 2020_hoogtestatistieken_gebouwen. Een workaround is om de tabel in de GeoPackage eerst te hernoemen met ogr2ogr (want ogr2ogr is niet zo kritisch) om dan je SQL query uit te voeren, zie het volgende Bash script:
curl -L https://download.pdok.nl/kadaster/basisvoorziening-3d/v1_0/2020/hoogtestatistieken/40bz1_2020_hoogtestatistieken_gebouwen.zip -o 40bz1_2020_hoogtestatistieken_gebouwen.zip \
&& unzip -o 40bz1_2020_hoogtestatistieken_gebouwen.zip \
&& ogr2ogr -f GPKG 40bz1_2020_hoogtestatistieken_gebouwen_renamed.gpkg 40bz1_2020_hoogtestatistieken_gebouwen.gpkg 2020_hoogtestatistieken_gebouwen -nln hoogtestatistieken_gebouwen_2020 \
&& ogrinfo 40bz1_2020_hoogtestatistieken_gebouwen_renamed.gpkg -sql "select * from hoogtestatistieken_gebouwen_2020 where Within(MakePoint(190922,443528),geom)"
Voor het uitvoeren van een SQL query met de ogr library in Python zijn genoeg voorbeelden te vinden op gis.stackexchange.
Wat betreft de niet geldige GeoPackage in het 3D basis bestand, ik ga hier intern melding van maken.
Voor andere opties voor het opvragen van hoogte zie dit draadje met bijv. :
- suggestie van @sbjager om hoogte informatie op te vragen met het AHN3 WMS GetFeatureInfo request
- mijn QGIS processing tool geschreven in Python voor het ophalen van hoogte informatie in QGIS
[edit]: Laatste versie van de processing tool staat in het pdokservicesplugin repo, ik ben bezig samen met @raymondnijssen en @rduivenvoorde met een nieuwe versie (met een PDOK processing toolbox).
Wat bedoel je hier precies met “locaties”? Je hebt de geometrie, da’s de locatie, dus ik begrijp niet helemaal wat je hiermee bedoelt.
Als het puur om een locatie in x en y gaat, is de WCS waarschijnlijk nog wel het makkelijkst, maar ik krijg het gevoel dat je met locatie iets anders bedoeld… een adres misschien?
1 like
Het gaat me inderdaad puur om de hoogte bepalen op een locatie (x, y). Ik zou uiteindelijk graag een script schrijven waar ik een x en een y coördinaat ingeef en dan als output een waarde krijg voor de hoogte op die locatie. Klinkt alsof het gebruik van de WCS een goede optie zou kunnen zijn.
Bedankt voor dit uitvoerige antwoord. Ik ga de verschillende opties uitzoeken. Klopt het dat ik de QGIS processing tool kan runnen in mijn IDE of moet ik deze runnen in de QGIS python console? En wat betreft WCS, klopt het dat ik een GetCoverage request kan doen met format image/tiff en dat ik op deze manier een GeoTIFF bestand terug krijg met hoogte waardes voor het gehele opgevraagde raster? Dit klinkt voor mij dan namelijk als de makkelijkste optie.
De processing tool die @antonbakker bedoelt zit als algoritme in de Processing Toolbox van QGIS. Deze is echter nog niet beschikbaar in de huidige PDOK Services plugin maar wordt hopelijk snel gereleased.
Ziet er dan zo uit:
Je kan de processing tool ook al handmatig installeren in QGIS; zie de instructies hier in de README.
Dat klopt, daarbij moet je wel rekening houden met de volgende twee zaken (zie ook deze reactie):
- de WCS service kan geen raster serveren van 1x1 pixel (bekende bug)
- als je een WCS 2.0.1 request doet krijg je een multipart response terug waar een geotiff in zit, wil je direct een geotiff response gebruik dan WCS versie 1.0.0 (voorbeeld 1.0.0 GetCoverage request). Het capabilities document geeft aan welke WCS versie’s ondersteund worden:
Dat gezegd hebbende, als je de hoogte van een gebouw wilt bepalen met de AHN3, zal je nog wel de terrein hoogte van de gevonden waarde af moeten trekken op die positie. Aangezien de AHN de absolute hoogte bevat t.ov. NAP. Zie ook de uitleg over het verschil tussen de DTM en DSM coverages in de AHN3 op de download pagina.
Persoonlijk zou ik de PDOK 3D hoogtestatistieken gebruiken om gebouw hoogte op te halen. Dan hoef je niet eerst te bepalen hoe hoog het maaiveld is rond het gebouw om de hoogte van het gebouw te bepalen.
Als ik de volgende code gebruik om een request te doen krijg ik de daaropvolgende response. Deze lijkt geen hoogte informatie te bevatten en dus niet correct te zijn. Kan je me vertellen wat ik fout doe?
Code:
import requests
response = requests.get(‘https://geodata.nationaalgeoregister.nl/ahn3/wcs?service=WCS&Request=GetCoverage&version=2.0.1&CoverageId=ahn3_05m_dsm&format=image/tiff&subset=x(196396.75,196397.75)&subset=y(450030.75,450032.75)’)
print(response.text)
Response:
–wcs
Content-Type: image/tiff
Content-Description: coverage data
Content-Transfer-Encoding: binary
Content-ID: coverage/out.tif
Content-Disposition: INLINE; filename=out.tif
II* �e(S�&��>��$n��a��� ���3.40282346638528898e+38�?�?f�aA�weA�����#��a@q)#a��+婁@Ӽ�I@J+�}@I�� ڿ�0(�hr�?w-!�������&S@Amersfoort / RD New|Amersfoort|1��BA�B���B��B
–wcs
Content-Type: application/octet-stream
Content-Description: coverage data
Content-Transfer-Encoding: binary
Content-ID: coverage/out.tif.aux.xml
Content-Disposition: INLINE; filename=out.tif.aux.xml
–wcs–
Dat is een multi-part antwoord, waarvan het 8-pixel grote tifje tussen de lege regel en de --wcs markering staat (1e deel, het begin zijn de mime-type definities). De overige data (2e deel, xml die je of niet meegekopieerd hebt, of die opgevreten is door het forum - want omgeven met xml-tags) is metadata. Je zul;t het dus moeten splitsen, en dan het tif-gedeelte binair opslaan.
In het xml-gedeelte van het antwoord staat de NoDataValue.
Edit: zie ook deze post
Je kunt dus het beste voor versie 1.0.0 gaan, zoals Anton ook al zegt, dan krijg je alleen een tifje terug.
1 like
hoe werkt dit aanpassen van de versie (Het lukt me niet het voorbeeld 1.0.0 GetCoverage request van Anton te openen)? Als ik in mn request version=2.0.1 vervang door version=1.0.0 krijg ik de volgende foutmelding. Ik heb in de foutmelding > vervangen door _ zodat het niet gefilterd wordt door het forum.
<?xml version='1.0' encoding="UTF-8" ?_ <ServiceExceptionReport version="1.2.0" xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/ogc https://geodata.nationaalgeoregister.nl/wcs/1.0.0/OGC-exception.xsd"_ <ServiceException code="MissingParameterValue" locator="crs"_msWCSGetCoverage(): WCS server error. Required parameter CRS was not supplied. </ServiceException_ </ServiceExceptionReport_Anton geeft een voorbeeld in zijn post.
voeg
crs=EPSG:28992&response_crs=EPSG:28992
toe…
PS. Als je html of xml code wilt toevoegen hier, dan kun je het plakken in de editor, selecteren, en als Vooraf opgemaakte tekst markeren (icoontje met </>). OP die manier kun je de xml hier wel zichtbaar krijgen.
Bedankt voor alle hulp op mijn vele vragen, maar met toevoeging krijg ik nog steeds een error.
https://geodata.nationaalgeoregister.nl/ahn3/wcs?service=WCS&Request=GetCoverage&version=1.0.0&CoverageId=ahn3_05m_dsm&format=image/tiff&subset=x(196396.75,196397.75)&subset=y(450030.75,450032.75)&crs=EPSG:28992&response_crs=EPSG:28992
<?xml version='1.0' encoding="UTF-8" ?>
<ServiceExceptionReport version="1.2.0"
xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/ogc https://geodata.nationaalgeoregister.nl/wcs/1.0.0/OGC-exception.xsd">
<ServiceException code="MissingParameterValue" locator="bbox/time">msWCSGetCoverage(): WCS server error. One of BBOX or TIME is required
</ServiceException>
</ServiceExceptionReport>
Is er misschien ergens documentatie te vinden over hoe een wcs request opgemaakt moet worden en welke dingen vereist zijn voor een verzoek?
Michiel gaf je dit draadje, waar heel veel info is te vinden in specifiek deze post: Starten met AHN WCS service - #2 door daneng - AHN - Geoforum
Onder andere een link naar de specs.
Anton gaf een voorbeeld van een 1.0.0 WCS Request:
https://geodata.nationaalgeoregister.nl/ahn3/wcs/v1_0?service=WCS&version=1.0.0&request=GetCoverage&coverage=ahn3_05m_dsm&crs=EPSG:28992&response_crs=EPSG:28992&bbox=214355,486976,214356,486977&format=image/tiff&resx=0.5&resy=0.5
Lijkt me dat je er daarmee toch wel moet komen? Er is een verschil in andere query parameters, alleen de toevoeging die ik aangaf is niet voldoende. Had ik wat beter moeten aangeven, sorry.
2 likes
ah vandaar. Heel erg bedankt. Hiermee kan ik aan de slag