Toch maar weer eens aan het spelen geweest met pc6 vlakken afleiden.
Ik neem nu spoorwegen/spoorwegemplacementen mee door deze uit de buurt polygonen te snijden (resultaat in de buurt_diff tabel). De ST_Intersection stap haalt dan ook die uit de postcodevlakken. Door de postcodevlakken op te splitsen en alleen die te recombineren die ook echt een adreslocatie hebben vallen ook de stukken eruit die ‘door/afgesneden’ worden door bv een spoorlijn.
Ik heb geprobeerd dit ook te doen met waterwegen en (hoofd/snel)wegen, maar daar werd het niet beter van naar mijn mening. Bij water zorgt het er bijvoorbeeld voor dat woonboten wegvallen. Dus voor nu maar even beperkt tot spoorwegen.
Wellicht is er nog ‘winst’ te halen door het type onderliggend terrein mee te nemen en wat weg te halen als er op “ander” terrein in zelfde vlak geen locatie is? Genoeg ideeën om verder te spelen 
.
Om de vlakken wat strakker te krijgen heb ik ook nog een ST_CoverageSimplify stap toegevoegd.
Informatie over spoorwegen haal ik uit de top10nl geopackage (terrein met typelandgebruik=‘spoorbaanlichaam’).
Volle SQL (die werkt op bag.duckdb, tabel buurt_diff zit daar ook in) is onderaan te vinden. Op mijn PC (Ryzen 5 2600 met 16GB) duurt alles zo’n beetje 3 minuten.
Resultaat is weer hier te vinden in pc6.gpkg
DuckDB SQL voor genereren PC6
-- 13s
-- 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,
st_makepoint(rd_x,rd_y)::GEOMETRY::GEOMETRY('EPSG:28992') AS loc
FROM addresses
WHERE postcode is not NULL
GROUP BY ALL;
-- 11s
-- create voronoi diagram and extract individual polygons
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;
-- 12s
-- Ensure everything is valid, otherwise next step may fail
update pc_vor set geom=st_makevalid(geom);
-- 26s
-- clip tiles to buurt_diff since DuckDB does not
-- support providing outer geom for ST_VoronoiDiagram
-- buurt_diff contains neighborhoods minus
-- any terrain in use as a railway.
create or replace table pc_vor_clipped as select v.bu_code, st_intersection(v.geom, bu.geom) as geom from pc_vor v left join buurt_diff bu on v.bu_code=bu.bu_code;
-- 5s
-- split multipolygons (bisected by step above)
-- so they fall out when combining with locations
-- in the next step
create or replace table pc_vor_clipped_split as select * replace(unnest(st_dump(geom)).geom as geom) from pc_vor_clipped;
-- 15s
-- Combine tiles with locations
create or replace table pc_com as select v.bu_code, v.geom, l.loc, l.min_pc,l.mode_pc as mode_pc, l.adr_cnt from pc_vor_clipped_split v inner join pc_loc l on ST_Within(l.loc, v.geom);
-- 1m15s
-- combine all tiles with the same postcode
-- and count the number of separate geometries
CREATE OR REPLACE TABLE pc6_raw AS
SELECT
bu_code,
mode_pc as pc6,
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;
-------------------------------
---- straigthen edges steps
-------------------------------
-- 2s
-- Split pc6 areas into separate parts (if they are)
create or replace table pc6_diss as select bu_code, pc6, unnest(st_dump(geom)).geom as geom from pc6_raw;
-- 6s
-- now simplify the shapes within coverages of a buurt/neighborhood
create or replace table pc6_smplfd as select bu_code, unnest(ST_Dump(ST_CollectionExtract(ST_CoverageSimplify(list(geom),10,false)))).geom as geom from pc6_diss group by all;
-- 11s
-- recombine with pc6
-- take these from pc_loc to find most appearing pc6
-- (borders may have shifted in the simplify step)
create or replace table pc6_smplfd_rec as select s.bu_code, s.geom, mode(l.mode_pc) as pc6 from pc6_smplfd s inner join pc_loc l on st_contains(s.geom, l.loc) group by all;
-- 5s
-- recombine al pc6 geoms again
CREATE OR REPLACE TABLE pc6 AS
SELECT
-- s.bu_code,
s.pc6,
p.address_count,
p.parts,
ST_Union_Agg(s.geom) AS geom
FROM pc6_smplfd_rec s
left join pc6_raw p on s.pc6=p.pc6
GROUP BY ALL;
.
.