2022 / Day 18: Blue
Today I’ll be looking at calculating shared paths of polygons with deriving left/right side information for them. This is useful for example for generating lines from admin area polygons where every border line is present only once. Or for example simplifying polygons without losing topological integrity.
I wrote about this same subject some time ago but using a different approach to do it. Time has passed and i’ve developed (hopefully) my skills and understanding, so I’m ready to try out st_sharedpaths.
My original approach had more aggregation in it so so hopefully this way it
will be more performant. As a test dataset I’ll generate a set of random points
with st_generatepoints
as in many previous times during this years #30DayMapChallenge, and cluster them
using st_clusterkmeans,
followed by doing a
st_voronoipolygons
to create a space evenly filled in with polygons (see CTE called voros).
As noted before - my original approach to solving this was to aggregate all
polygon boundaries and then st_linemerge
them. Which can become a rather tedious task if there’s a lot of data.
A better way around is with st_sharedpaths
which works with a pair of geometries and returns a GeometryCollection
of two elements:
- the first one is lines going in the same direction,
- the second is lines going in the opposite directions
In order to rule out duplicates I’ll condense the pairing geometries into
and array ordered by the value of their respective identifiers (cl in
this case) which allows me to do a
select
distinct on (<ordered_identifier_array>)
<ordered_identifier_array>, <ordered_geoms_array>
from
...
Another way for getting unique combinations here would be to use
select
*
from (
select
row_number() over(partition by <ordered_identifier_array>) as rn
<ordered_identifier_array>, <ordered_geoms_array>
from
...
) x
where
x.rn = 1
as I used for 2022 / Day 6: Network. The result of
this is encapsulated in CTE pairs.
For attaching L/R side info to the linestring I previously used an approach of doing a minuscule L/R st_offsetcurve followed by st_lineinterpolatepoint to the midway of the offset linestring and then st_within with the midway point to discover what’s on either side of the line. Now with st_sharedpaths it’s a bit simpler because:
- out of the pair of geometries we know which one is first and which one second
- st_sharedpaths returns linestrings in the direction of the first geometry.
So unless the polygon interiors share space / overlap (meaning their shared paths of
boundaries are moving in the same direction) we can safely assume that
the first of the geometry pair is always to the right, and the second one to
the left for clockwise (st_ispolygoncw)
and the other way around for counter-clockwise
(st_ispolygonccw) polygons.
Which is also the reason why
st_forcepolygoncw is used after
the voronoi cell unioning in CTE voros. Although it doesn’t do anything in
this case - our polygons will be clockwise by default.
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
),
voros as (
select cl, ints as geom
from
bounds, (
select
c.cl,
st_forcepolygoncw(st_union(d.geom)) as geom
from (
select
st_voronoipolygons(
st_collect(geom)
) as geom
from
clusters
) v
join lateral
st_dump(v.geom) d on true
join lateral (
select c.cl
from clusters c
where st_within(c.geom, d.geom)
limit 1
) c on true
group by c.cl
) b
join lateral
st_intersection(
b.geom,
bounds.geom
) ints on true
),
pairs as (
select
distinct on (cls)
cls, geoms
from (
select
a.cl,
case
when a.cl < b.cl then
array[a.cl, b.cl]
else
array[b.cl, a.cl]
end as cls,
case
when a.cl < b.cl then
array[a.geom, b.geom]
else
array[b.geom, a.geom]
end as as geoms
from
voros a,
voros b
where
st_intersects(a.geom, b.geom) and
a.cl != b.cl
) d
),
shared_paths as (
select
row_number() over()::int as oid, pairs.cls,
st_linemerge(st_collect(d.geom)) as geom
from pairs
join lateral
st_sharedpaths(
st_boundary(geoms[1]),
st_boundary(geoms[2])
) sp on true
join lateral
st_dump(sp) d on true
group by pairs.cls
)
select
sp.oid,
sp.cls[1] as right_cl, sp.cls[2] as left_cl,
sp.geom
from
share_paths sp
;