code

PDOK BGT OGC API Feature collections

Python script om inzicht te krijgen welke velden en veld type via de API terug kunnen worden gegeven en in welk veld de waarde 'brievenbus' terugkomt.

###
import requests

# URL om de data op te halen

url = "https://api.pdok.nl/lv/bgt/ogc/v1/collections/straatmeubilair/items?f=json&limit=1000"

# Maak een verzoek naar de API

response = requests.get(url)

# Controleer of het verzoek succesvol was

if response.status_code == 200:

    data = response.json()

    # Haal de eerste feature op om de structuur te analyseren

    if data['features']:

        first_feature = data['features'][0]

        properties = first_feature['properties']

        # Print veldnamen en datatypes

        print("Veldnamen en datatypes van de BGT straatmeubilair features:")

        for key, value in properties.items():

            print(f"{key}: {type(value).__name__}")

        # Bepaal in welk veld de waarde 'brievenbus' wordt opgeslagen

        for feature in data['features']:

            if 'brievenbus' in feature['properties'].values():

                for key, value in feature['properties'].items():

                    if value == 'brievenbus':

                        print(f"De waarde 'brievenbus' wordt opgeslagen in het veld: {key}")

                        break

                break

    else:

        print("Geen features gevonden.")

else:

    print(f"Error fetching data: {response.status_code}")

    print(f"Response content: {response.text}")

###

PDOK BGT OGC API Feature selectie

Python script om alle brievenbussen in de gemeente Den Haag uit de BGT data te halen via de PDOK OGC API Features en vervolgens weg te schrijven in een feature class in ARCGIS Pro 3.2 om deze vervolgens op de kaart te tonen.

###
import arcpy

import requests

import json

# Set up the environment

arcpy.env.overwriteOutput = True

workspace = r"C:\Users\rv\Documents\ArcGIS\Projects\PDOK\PDOK.gdb"

out_fc = "brievenbussen_BGT"

# Define spatial references

sr_wgs84 = arcpy.SpatialReference(4326)  # WGS84

sr_amersfoort = arcpy.SpatialReference(4289)  # Amersfoort

sr_rd_new = arcpy.SpatialReference(28992)  # RD New

# API endpoint

url = "https://api.pdok.nl/lv/bgt/ogc/v1/collections/straatmeubilair/items"

# Parameters for Den Haag bounding box (approximate)

params = {

    "bbox": "4.21,52.01,4.41,52.13",

    "limit": 1000,

    "f": "json"

}

# Function to get features from a single page

def get_features(url, params):

    response = requests.get(url, params=params)

    if response.status_code == 200:

        data = response.json()

        return data['features'], data.get('links', [])

    else:

        print(f"Error: {response.status_code}")

        return [], []

# Fetch all features using cursor pagination

all_features = []

while True:

    features, links = get_features(url, params)

    all_features.extend(features)

    next_link = next((link for link in links if link['rel'] == 'next'), None)

    if next_link:

        url = next_link['href']

        params = {}  # Reset params as the new URL includes them

    else:

        break

print(f"Total features retrieved: {len(all_features)}")

# Filter features to only include 'brievenbus'

filtered_features = [feature for feature in all_features if feature['properties'].get('plus_type') == 'brievenbus']

print(f"Total 'brievenbus' features: {len(filtered_features)}")

# Create a new feature class

arcpy.management.CreateFeatureclass(workspace, out_fc, "POINT", spatial_reference=sr_rd_new)

# Add fields to the feature class

arcpy.management.AddFields(out_fc, [

    ["identificatie", "TEXT"],

    ["bgtType", "TEXT"],

    ["plus_type", "TEXT"],

    ["relatieveHoogteligging", "SHORT"],

    ["status", "TEXT"],

    ["x_coord", "DOUBLE"],

    ["y_coord", "DOUBLE"]

])

# Insert features into the feature class

with arcpy.da.InsertCursor(out_fc, ["SHAPE@", "identificatie", "bgtType", "plus_type", "relatieveHoogteligging", "status", "x_coord", "y_coord"]) as cursor:

    for feature in filtered_features:

        properties = feature['properties']

        geometry = feature['geometry']

        if geometry['type'] == 'Point':

            x, y = geometry['coordinates']

            # Create point geometry in WGS84

            point_wgs84 = arcpy.PointGeometry(arcpy.Point(x, y), sr_wgs84)

            # Project to Amersfoort

            point_amersfoort = point_wgs84.projectAs(sr_amersfoort)

            # Project to RD New

            point_rd_new = point_amersfoort.projectAs(sr_rd_new)

            # Get x and y coordinates in RD New

            x_rd, y_rd = point_rd_new.firstPoint.X, point_rd_new.firstPoint.Y

            # Extract plus_type, defaulting to None if not present

            plus_type = properties.get('plus_type')

            cursor.insertRow([

                point_rd_new,

                properties.get('identificatie'),

                properties.get('bgtType'),

                plus_type,

                properties.get('relatieveHoogteligging'),

                properties.get('status'),

                x_rd,

                y_rd

            ])

print(f"Feature class '{out_fc}' created with {len(filtered_features)} features.")
###

Overpass API selectie

Python script om alle brievenbussen in de gemeente Den Haag uit de OpenStretMap data te halen via de Overpass API en vervolgens weg te schrijven in een feature class in ARCGIS Pro 3.2 om deze vervolgens op de kaart te tonen.

###
import arcpy

import requests

import json

# Set up the environment

arcpy.env.overwriteOutput = True

workspace = r"C:\Users\renev\Documents\ArcGIS\Projects\PDOK\PDOK.gdb"

out_fc = "brievenbusen_OSM"

# Define spatial references

sr_wgs84 = arcpy.SpatialReference(4326)  # WGS84

sr_rd_new = arcpy.SpatialReference(28992)  # RD New

# Overpass API endpoint

overpass_url = "http://overpass-api.de/api/interpreter"

# Query to get all mailboxes in Den Haag

overpass_query = """

[out:json][timeout:25];

(

  node["amenity"="post_box"](52.01,4.21,52.13,4.41);

);

out body;

>;

out skel qt;

"""

# Send request to Overpass API

response = requests.get(overpass_url, params={'data': overpass_query})

data = response.json()

# Create a new feature class

arcpy.management.CreateFeatureclass(workspace, out_fc, "POINT", spatial_reference=sr_rd_new)

# Add fields to the feature class

arcpy.management.AddFields(out_fc, [

    ["osm_id", "TEXT"],

    ["ref", "TEXT"],

    ["operator", "TEXT"]

])

# Insert features into the feature class

with arcpy.da.InsertCursor(out_fc, ["SHAPE@", "osm_id", "ref", "operator"]) as cursor:

    for element in data['elements']:

        if element['type'] == 'node':

            x, y = element['lon'], element['lat']

            # Create point geometry in WGS84

            point_wgs84 = arcpy.PointGeometry(arcpy.Point(x, y), sr_wgs84)

            # Project to RD New

            point_rd_new = point_wgs84.projectAs(sr_rd_new)

            # Extract attributes

            osm_id = str(element['id'])

            ref = element['tags'].get('ref', '')

            operator = element['tags'].get('operator', '')

            cursor.insertRow([

                point_rd_new,

                osm_id,

                ref,

                operator

            ])

print(f"Feature class '{out_fc}' created with {len(data['elements'])} mailboxes.")

###