Vertalen AHN3 tegels naar Rijksdriehoek coordinaten

Datasets AHN
#1

Hallo,

Ik ben benieuwd hoe ik de AHN3 tegelnamen kan converten naar X en Y Rijksdriehoek coordinaten (bijvoorbeeld de min x en min y RD coordinaten van een tegel).

Na het downloaden van een aantal AHN tegels heb ik de BBOX coordinaten kunnen berekenen met de volgende code:

from geotiff import GeoTiff
# Get the bounding box
bbox = geotiff.tif_bBox
x_min, y_max = bbox[0]
x_max, y_min = bbox[1]

En de output in csv:

tile_name,x_min,y_max,x_max,y_min
25BZ1,110000.0,493750.0,115000.0,487500.0
25BZ2,115000.0,493750.0,120000.0,487500.0
25EZ1,120000.0,493750.0,125000.0,487500.0
25EZ2,125000.0,493750.0,130000.0,487500.0
25FZ1,130000.0,493750.0,135000.0,487500.0
25DN1,110000.0,487500.0,115000.0,481250.0
25DN2,115000.0,487500.0,120000.0,481250.0
25GN1,120000.0,487500.0,125000.0,481250.0
25GN2,125000.0,487500.0,130000.0,481250.0
25HN1,130000.0,487500.0,135000.0,481250.0
25DZ1,110000.0,481250.0,115000.0,475000.0
25DZ2,115000.0,481250.0,120000.0,475000.0
25GZ1,120000.0,481250.0,125000.0,475000.0
25GZ2,125000.0,481250.0,130000.0,475000.0
25HZ1,130000.0,481250.0,135000.0,475000.0

Is het ook mogelijk om deze RD coordinaten te berekenen met alleen de naam van een AHN3 tegel?

#2

Je kunt misschien beter de Bladindex WFS inlezen in QGIS en van daaruit een export draaien:

https://geodata.nationaalgeoregister.nl/ahn3/wfs?request=GetCapabilities

image
image819×924 27.9 KB

Aan de bladnaam kun je niks afleiden. Zo zijn er wat aanvullende bladen bij Limburg dacht ik en de tweede maasvlakte bijgekomen zonder dat de overige bladen hierop hernummerd zijn.

1 like
#3

Ter aanvulling op @Anton, hier is nog een voorbeeld Python script om de WFS GetFeature reponse om te zetten in CSV output:

#!/usr/bin/env python3
from urllib.request import urlopen
import json
import sys
import csv

def get_coordinates_features(features):
    def map_fun(feature):
        minx = feature["geometry"]["coordinates"][0][0][0][0] # multipolygon causing the nesting
        miny = feature["geometry"]["coordinates"][0][0][0][1]
        maxx = feature["geometry"]["coordinates"][0][0][2][0]
        maxy = feature["geometry"]["coordinates"][0][0][2][1]
        bladnr = feature["properties"]["bladnr"]
        return bladnr,minx,miny,maxx,maxy
    return [map_fun(feature) for feature in features]
    
def main():
    url = "https://geodata.nationaalgeoregister.nl/ahn3/wfs?request=GetFeature&typename=ahn3:ahn3_bladindex&outputFormat=json"
    with urlopen(url) as response:
        json_string = response.read()
    data = json.loads(json_string)
    coordinates = get_coordinates_features(data["features"])
    csv_writer = csv.writer(sys.stdout)
    csv_writer.writerow(["bladnr", "minx","miny","maxx","maxy"])
    for coord in coordinates:
        csv_writer.writerow(coord)
    
if __name__ == "__main__":
    main()
1 like
#4

Dit gaat over de bladindeling in plaats van tegels. Omdat de nummering van de bladen aangepast is op de vorm van Nederland kan je de begrenzingen niet uitrekenen zoals je met tegels volgens de Nederlandse tiling-richtlijn (zoals https://geodata.nationaalgeoregister.nl/ahn3/wms?request=GetCapabilities&service=wms) wel kan. Zie voor die formule: Vertalen tiling (XYZ tiles) naar Rijksdriehoek coordinaten - PDOK Kaart - Geoforum