2022 / Day 02: Lines

Moving on with generative SQL. This exercise attempts to create a set of linestrings which would connect nodes on all four sides of a square box ( in this instance as well as yesterday i was using the bbox of epsg:3301 based tile gridset) so that all nodes would be connected with every other node except the ones on the same side.

This time around there’s a lot of lateral joins, and in real life most of these CTE queries should/could be written to temporary tables in order to enforce indexes etc. But staying with the spirit of it has to run as a single statement, well… yes.

Some (but not all) PostGIS functions that are doing the heavy lifting here:

st_makeline to construct a line based on a bunch of consecutive points. Although there are variants of two-point-makeline and aggregate-makeline, over the years i’ve really grown accustomed to st_makeline(array[a.geom, b.geom]) and st_makeline(array_agg(a.geom order by a.i)). Maybe because using arrays seems to make everything fancier?
st_dumppoints to extract all nodes of a geometry. NB! pay attention to the composition of the returned geometry_dump.path composition. This will vary depending the type of geometry fed into the function
st_segmentize to create more nodes on a linestring (or polygon boundary for that matter) so that no segment (linestring bit between two consecutive nodes) is longer than the input distance.

The basic idea behind this query is to create a square box, dump the corner nodes with st_dumppoints and create segments (of two consecutive nodes of the ring).

NB! As of PostGIS 3.2 we can use st_dumpsegments instead of

select st_makeline(array[pt.geom, lead(pt.geom) over (order by pt.path)]
from st_dumppoints(geom) pt)

to achieve the same thing.

Now for every returned segment, lets add extra points every 10000 meters with st_segmentize, and st_dumppoints to extract all points of the box.

Select unique nodes by st_union followed by st_dump to get node geometries and identifiers.

For every box point find the closest node and it’s identifier with a lateral join and st_dwithin selecting the closest with

order by nodes.geom <-> box_points.geom
limit 1

And finally all that’s left to do is to find all possible combinations of nodes, and st_makeline out of these. I’m using ordering by node_id value so I can select only unique lines. Meaning if there’s a line going from A->B then I don’t want/need the reversed B->A in the output.

And done…

output image

with
    /* minmax defines the corners where we create a set of random points. */
    minmax as (
        select
            st_point(40500.000000,5993000.000000,3301) as ll,
            st_point(1064500.000000,7017000.000000,3301) as ur
    ),
    /* bounds gives the envelope for them  */
    bounds as (
        select
            1 as id,
            st_envelope(
                st_collect(
                    array[ll, ur]
                )
            ) as geom
        from minmax
    ),
    /* extract envelope sides (4 pcs) */
    box_lines as (
        select
            bounds.id,
            pt.path[2] as p,
            st_makeline(
                array[
                    pt.geom,
                    lead(pt.geom) over (
                        partition by bounds.id, pt.path[1]
                        order by pt.path[2]
                    )
                ]
            )::geometry(linestring, 3301) as geom,
            min(x) over (partition by bounds.id, pt.path[1]) as minx,
            min(y) over (partition by bounds.id, pt.path[1]) as miny,
            max(x) over (partition by bounds.id, pt.path[1]) as maxx,
            max(y) over (partition by bounds.id, pt.path[1]) as maxy
        from
            bounds
                join lateral
                    st_dumppoints(bounds.geom) pt on true
                join lateral
                    st_x(pt.geom) as x on true
                join lateral
                    st_y(pt.geom) as y on true
    ),
    /* generate points - segmentize the box side line by 50 km */
    box_pts as (
        select
            id as box_id, p as box_side, minx, miny, maxx, maxy, x, y,
            pts.geom, pts.path[1] as ord,
            max(pts.path[1]) over (partition by id, p) as max_ord_on_side
        from
            box_lines
                join lateral
                    st_segmentize(geom, 50000) seg on true
                join lateral
                    st_dumppoints(seg) pts on true
                join lateral
                    st_x(pts.geom) x on true
                join lateral
                    st_y(pts.geom) y on true
        where
            p < 5
    ),
    /* generate "singular" nodes and unique node_ids */
    nodes as (
        select
            path[1] as node_id, geom as geom
        from (
            select
                (st_dump(st_union(box_pts.geom))).*
            from
                box_pts
        ) n
    ),
    /* merge box_pts with nodes - essentially assigns node_ids to box_pts */
    pts as (
        select
            p.box_id, p.box_side, p.geom, p.ord
            p.x, p.y, p.minx, p.miny, p.maxx, p.maxy,
            n.node_id
        from
            box_pts p
                join lateral (
                    select node_id
                    from nodes
                    where st_dwithin(nodes.geom, p.geom, 1)
                    order by nodes.geom <-> p.geom
                    limit 1
                ) n on true
    )
/* ... AND create the linework connecting every node to every other node
   unless the nodes are on the same side of the box.
*/
select
    row_number() over()::int as oid, box_id, node_ids, geom, st_length(geom)
from (
    select
        /* get distinct pairs regardless on direction by */
        --distinct on (node_ids)
        /* OR INSTEAD: */
        row_number() over (partition by box_id, node_ids) as x,
        box_id, node_ids, geom
    from (
        select
            a.box_id,
            case
                when a.node_id < b.node_id then st_makeline(array[a.geom, b.geom])
                else st_makeline(array[b.geom, a.geom])
            end as geom,
            case
                when a.node_id < b.node_id then array[a.node_id, b.node_id]
                else array[b.node_id, a.node_id]
            end as node_ids
        from
            pts a,
            pts b
        where
            a.box_id = b.box_id and
            a.box_side != b.box_side and (
                (
                    not (a.y = b.y or a.x = b.x ) or
                    (a.y = b.y and a.y != all(array[a.miny, a.maxy])) or
                    (a.x = b.x and a.x != all(array[a.minx, a.maxx]))
                ) or  (
                    a.ord = 1 and
                    b.ord = 1
                )
            )
    ) v
) b
where
    x = 1
;