Heatmaps - van postcode naar GPS coördinaten

De overheidsorganisatie PDOK - Publieke Dienstverlening Op de Kaart - houdt een registratie bij van alle objecten in Nederland. Dit is het zogeheten 'BAG'-register, Basisregistratie Adressen Gebouwen. Hierin staan de adressen, soort objecten en coordinaten van alle gebouwen in Nederland.

Dit BAG-register is te downloaden in de vorm van een SQLite-database van ~7gb. Voor de coördinaten worden extensies gebruikt, daarom is het handiger spatialite te gebruiken om het bestand te openen. Zo kan je direct bepaalde BLOB-velden inspecteren.

De locaties zijn genoteerd in zogeheten rijksdriehoek-coördinaten die kunnen worden opgezet naar GPS-coördinaten. Gelukkig is al iemand me voor geweest hier een library voor te maken, te vinden op github: https://github.com/djvanderlaan/rijksdriehoek/blob/master/Python/rijksdriehoek.py

Het gave is dat je vervolgens mee heatmaps kan maken van alle postcode/huisnummer bestanden!

#!/usr/bin/python3


# BAG (EPSG:28992) Geopackage @ https://service.pdok.nl/lv/bag/atom/bag.xml

# spatialite gui & @ https://www.gaia-gis.it/gaia-sins/ / https://www.gaia-gis.it/gaia-sins/windows-bin-amd64-latest/

# rijksdriehoek.py @ https://github.com/djvanderlaan/rijksdriehoek/blob/master/Python/rijksdriehoek.py


import sqlite3
from rijksdriehoek import rd_to_wgs


def dict_factory(cursor, row):
    fields = [column[0] for column in cursor.description]
    return {key: value for key, value in zip(fields, row)}


# connect to db
con = sqlite3.connect("C:/tmp/postcodes/bag-light.gpkg")
con.enable_load_extension(True)

con.execute('SELECT load_extension("mod_spatialite")')   
con.execute('SELECT InitSpatialMetaData(1);')  
con.row_factory = dict_factory



cur = con.cursor()

cur.execute("select (ST_MaxX(vo.geom)+ST_MinX(vo.geom))/2 rdx, (ST_MaxY(vo.geom)+ST_MinY(vo.geom))/2 rdy, vo.* from verblijfsobject vo")

for row in cur :
    #print( row )
    straat = row['openbare_ruimte_naam']
    huisnr = row['huisnummer']
    huisletter = row['huisletter']
    huistoevoeging = row['toevoeging']
    postcode = row['postcode']
    plaats = row['woonplaats_naam']
    
    latlon = rd_to_wgs(row['rdx'], row['rdy'])
    print( latlon, straat, huisnr, huisletter, huistoevoeging, postcode, plaats )