2022 / Day 09: Space
I wrote a little something about this same subject a few years ago in the context of expanding admin unit polygon boundaries out to the sea using coastline vertices. And for a nice homogeneous dataset it works really nicely. Now if you apply the same thing in a different kind of situation (e.g. filling in holes between cadastral parcels of varying sizes and most importantly with varying amounts of vertices on the sides facing but not quite touching each other) it might not be the most suitable solution.
In todays challenge I will try to tackle the problem of saw-teeth that arises from… too little points. But that said, adding more vertices to polygon boundary will only make the saw-teeth finer in resolution and will not deal with the problem.
Todays output image will be an animated gif showing the different stages (CTE results) that the query produces. The query itself (see below) is still a oneliner.
The query itself starts off again with defining the minmax corners, and
constructing a bounding geometry into bounds using
st_envelope. Then generating a
bunch of random points (st_generatepoints)
and clustering them with st_clusterkmeans
into clusters. Now in order to get concave polygons instead of convex ones
st_concavehull is used with
0.0 produce a hull of maximum concaveness as areas. For extra fun.
So far so good. Until now we’ve produced only the testing dataset. This is where we start the real processing.
First off, extract all area boundary vertices with
st_dumppoints,
and keeping in mind that the boundary linestring extracted with
st_boundary will have
first and last points in the same location so
filter last one out using the number of vertices per geometry calculated by
st_numpoints.
All of this is in the CTE called pts.
Secondly, we’ll also extract all segments of area boundaries with
st_dumpsegments into segs.
For some of them we’ll add extra vertices: for every point in pts find the
closest other segment (st_dwithin
and must not belong to the same area), locate the closest
point along the segment with
st_linelocatepoint (
which returns a fraction between 0 and 1) and then
st_lineinterpolatepoint
to get the geometry of the extra point for the other area (the one,
whose segment this is) which would balance out the effect of area A.node X
having too much influence on the neighboring area with no vertices in sight
nearby. Out of all of these possible extra_pts it makes sense to use only
the ones where a linestring constructed (
st_makeline) from a node on the
boundary of area A to an interpolated closest point location on the boundary
of area B does not st_crosses
the boundary of area A.
Then st_collect the points
(both pts and extra_pts) and then generate
st_voronoipolygons followed
by a st_dump to single part polygons
and discovering which area identifier (cluster_id in this case) the voronoi
cell belongs to by joining pts and extra_pts with
st_within.
Just to be on the safe side, lets do a
st_snaptogrid followed by
st_makevalid
and st_dump again so we can filter out
only POLYGON type geometries in the last step.
The last step consists of aggregating all voronoi cells with the same
cluster_id value using st_union and
then finding their intersection with the original bounds geometry
with st_intersection.
The result is still not perfect but i think the direction this is going is
producing way better results than previously. One thing to be thought about
and maybe elaborated a bit more is how to create the extra_pts on neighboring
areas boundaries. Currently only the closest one is picked. But as it shows up
this might be inadequate in a situation where multiple area corners are close
by. As the saying goes: more research needs to be done. :)
with
minmax as (
select 1 as oid,
st_point(
40500.000000,5993000.000000,3301
) as ll,
st_point(
1064500.000000,7017000.000000,3301
) as ur
),
bounds as (
select
st_envelope(
st_collect(
array[ll, ur]
)
) as geom
from minmax
),
clusters as (
select
pt.path[1] as oid, pt.geom,
st_clusterkmeans(
pt.geom,
30,
10000000
) over () cl
from
bounds
join lateral
st_generatepoints(
bounds.geom,
1000
) pts on true
join lateral
st_dump(
pts
) pt on true
),
areas as (
select
cl as cluster_id, count(1),
st_concavehull(st_collect(geom),0.0) as geom
from
clusters
group by
cl
having
count(1) > 5
),
pts as (
select
cluster_id, pts.path[1] as pnt_ord,
pts.geom
from
areas
join lateral
st_boundary(areas.geom) b on true
join lateral
st_dumppoints(b) pts on true
where
pts.path[1] < st_numpoints(b)
),
segs as (
select
areas.cluster_id, seg.path[1] as p,
seg.geom
from
areas
join lateral
st_dumpsegments(
areas.geom
) seg on true
),
extra_pts as (
select
c.cluster_id as closest_cluster,
c.closest_point, c.connector
from
pts a
join lateral (
select
f.cluster_id, closest_point,
st_makeline(
a.geom,
closest_point
) as connector
from
areas, (
select
segs.cluster_id,
segs.geom,
st_linelocatepoint(
segs.geom,
a.geom
) as fr
from
segs
where
st_dwithin(
a.geom,
segs.geom,
100000
) and
a.cluster_id != segs.cluster_id
order by
segs.geom <-> a.geom
limit 1
) f
join lateral
st_lineinterpolatepoint(
f.geom,
f.fr
) closest_point on true
where
f.fr > 0.0 and
f.fr < 1.0 and
st_crosses(
st_makeline(
a.geom,
closest_point
),
areas.geom
) = false and
a.cluster_id = areas.cluster_id
) c on true
),
voro as (
select
coalesce(cl_pts.cluster_id, x_pts.cluster_id) as cluster_id,
s.path[1] as p, s.geom as geom
from (
select
st_voronoipolygons(
st_collect(
d.geom
)
) as geom
from (
select geom from pts
union all
select closest_point from extra_pts
) d
) v
join lateral
st_dump(v.geom) d on true
left join lateral (
select
cluster_id
from
pts
where
st_within(pts.geom, d.geom)
limit 1
) cl_pts on true
left join lateral (
select
closest_cluster as cluster_id
from
extra_pts
where
st_within(extra_pts.closest_point, d.geom)
limit 1
) x_pts on true
join lateral
st_dump(
st_makevalid(
st_snaptogrid(
d.geom,
0.1
)
)
) s on true
)
/* And union voronoi cells by their area identifier,
and intersection with bounds*/
select
row_number() over()::int as oid,
cluster_id, d.geom
from
bounds, (
select
cluster_id, st_union(geom) as geom
from
voro
group by
cluster_id
) v
join lateral
st_intersection(bounds.geom, v.geom) i on true
join lateral
st_dump(i) d on true
where
geometrytype(d.geom) = 'POLYGON'
;