Wat DuckDB niet ondersteunt is bijvoorbeeld een boundary geometry bij het creëren van een Voronoi diagram, dus dat moet dan in een aparte stap worden geclipt met een ST_Intersection. Dat doe ik nu dus aan de hand van buurten. Daar waar een postcode zo’n grens overgaat lijkt dat niet veel uit te maken, die sluiten dan wel weer aan (wellicht wat hoekiger); maar er zijn sowieso ook al zat postcodegebieden die niet aansluiten.
Bij woonplaatsen mis ik echt een paar best grote vlakken, het zijn geen randgevallen. Daar nu naar kijkend denk ik dat dit een fout is in mijn code bij het inlezen van de BAG. Het lijken allemaal woonplaatsen met disjuncte oppervlakten (Elim, Blauwestad, Noordseschut, Leerbroek, Geldermalsen, maar ook Amsterdam), dat zal helpen bij het zoeken van de oorzaak 
. Toch nuttig zo’n compleet ongerelateerde zoektocht 
.
Totale tijd van alle stappen (op mijn meer dan 6 jaar oude Ryzen 5 2600) is onder de 2,5 minuut. Zie beneden voor de SQL met timings.
De beschrijving van @JBak (dank!) geeft richting en kan ik mooi mee verder knutselen, dat is toch het leukst.
DuckDB SQL voor genereren postcode6
-- 17s
-- Determine all unique locations
-- Assign most common postcode at overlaps
CREATE OR REPLACE TABLE pc_loc AS
SELECT
bu_code,
min(postcode) AS min_pc,
list_mode(list_sort(list(postcode))) as mode_pc,
list_sort(list(DISTINCT postcode)) AS all_pc,
count(*) AS adr_cnt,
loc AS loc
FROM addresses
WHERE postcode is not NULL
GROUP BY bu_code, loc;
-- 10s
-- create voronoi diagram and unnest tiles
create or replace table pc_vor as select
bu_code,
unnest(ST_Dump(ST_VoronoiDiagram(ST_Collect(list(loc))))).geom as geom
from pc_loc
group by bu_code;
-- 37s
-- clip tiles to buurt since DuckDB does not
-- support providing outer geom for ST_VoronoiDiagram
create or replace table pc_vor1 as select
v.bu_code,
st_intersection(v.geom, bu.geom) as geom
from pc_vor v left join buurt bu on v.bu_code=bu.bu_code;
-- 12s
-- Combine tiles with locations
-- appr. 1800 locations will have no geom!
-- (at the edge of gemeente, e.g. near coast/water)
create or replace table pc_com as select
v.bu_code,
v.geom,
l.loc,
l.min_pc,
l.mode_pc,
l.adr_cnt
from pc_loc l left join pc_vor1 v on ST_Within(l.loc, v.geom);
-- 1m5s
-- combine all tiles with the same postcode
-- and count the number of separate geometries
CREATE OR REPLACE TABLE pc6 AS
SELECT
mode_pc as postcode6,
sum(adr_cnt)::INTEGER AS address_count,
ST_NGeometries(ST_Union_Agg(geom)) AS parts,
ST_Union_Agg(geom) AS geom
FROM pc_com
GROUP BY all;
.