Hoogtewaarden AHN DTM_05

beste allemaal,

Ik ben bezig met een project van de TU Delft en ik krijg die ahn kaartbladen niet gelezen via python. ik moet de hoogtes van rotterdam hebben. Van kaartnr: 37GN2, 37EZ2, 37FZ1, 37HN1. Elke keer staat er dit: RasterioIOError: Read failed. See previous exception for details. In mijn code wordt ook kaartblad verwerkt, maar dat duurt heel lang, dus vroeg mij af of er een andere manier is om in python de hoogtes te weten per coordinaat. Ik hoor het graag.

dit is mijn code:

import rasterio
import rasterio.features
import rasterio.warp
import requests
from io import BytesIO

Dictionary met kaartbladen en bijbehorende URLs

kaartbladen = {
“M_37GZ1”: “https://service.pdok.nl/rws/ahn/atom/downloads/dtm_05m/M_37GZ1.tif”,
“M_37FZ1”: “https://service.pdok.nl/rws/ahn/atom/downloads/dtm_05m/M_37FZ1.tif”,
“M_37HN2”: “https://service.pdok.nl/rws/ahn/atom/downloads/dtm_05m/M_37HN2.tif”,
“M_37EN1”: “https://service.pdok.nl/rws/ahn/atom/downloads/dtm_05m/M_37EN1.tif
}

for naam, url in kaartbladen.items():
print(f"Verwerken van kaartblad: {naam}")

# Download het .tif bestand van de URL
print(f"Downloaden van {naam}...")
response = requests.get(url)
if response.status_code == 200:
    print(f"Succesvol gedownload {naam}, bezig met openen van het bestand...")
    with rasterio.open(BytesIO(response.content)) as dataset:
        print(f"Bestand geopend voor {naam}, bezig met het laden van masker...")
        mask = dataset.dataset_mask()

        # Check of het masker leeg is
        if mask.sum() == 0:
            print(f"Geen geldige data gevonden in kaartblad {naam}")
            continue

        print(f"Masker geladen voor {naam}, bezig met het ophalen van geometrieën...")

        # Haal geometrische vormen en waarden op uit de array
        count = 0
        for geom, val in rasterio.features.shapes(mask, transform=dataset.transform):
            count += 1
            if count % 100 == 0:  # Print om de 100 geometrieën
                print(f"{count} geometrieën verwerkt voor {naam}")

            # Transformeer de vormen naar EPSG:4326
            geom = rasterio.warp.transform_geom(dataset.crs, 'EPSG:4326', geom, precision=6)
            
            # Print GeoJSON-vormen naar stdout
            print(geom)
        print(f"Verwerking van {naam} voltooid, {count} geometrieën verwerkt.")
else:
    print(f"Fout bij het downloaden van {naam}: {response.status_code}")

Geen direct antwoord op je vraag maar ik denk dat dit heel mooi zou zijn. Heb je deze python code bekeken? Die haalt ruwe AHN scandata op.Ik zou AHN_05 in de .laz of .csv willen krijgen. .tif heeft de ontwerper niet veel aan omdat het plat 2D is.

HideBa/ahn_cli: Python CLI tool which helps you to download AHN point cloud data easily (github.com)

video’tje

Oo, en deze code is intressant

ahninpainter/utils/download.py at main · chenzhaiyu/ahninpainter (github.com)

Welke geometrie verwacht je aan te treffen (en om te zetten naar WGS84)?

Overweeg het probleem sterk te versimepelen door de data eerst te downloaden en daarna lokaal te gebruiken in plaats van alles in één keer te doen.