2022 / Day 12: Scale

Let’s run another comparison of PostGIS functions. This time simplification which is a an integral part of preparing data for maps on multiple scales.

We’ll look at:

For the data we’ll use a linestring in a shape of “S”, and generate_series a few copies of them into boxes of equal width and height: 400 m * 400 m. Along the y- axis of the plot we’ll have different simplification functions battling out, and along x different levels of tolerance.

To ease the understanding of what’s going on, every box has two numbers in the lower right corner. The first one (upper) is the number of nodes that are left to the geometry after simplification with the function. The lower one shows the average length of segments (derived with st_dumpsegements) from the simplification.

Thing to note here is that (as the name says) st_reduceprecision works a bit differently than the rest - its main job is to reduce the precision of coordinates, much like st_snaptogrid only that it will not return invalid geometries and will remove those completely that are below the tolerance.

output image

with
    col as (
        /* generate a column*/
        select
            y,
            st_translate(geom, 0, y*400) as geom,
            st_translate(box, 0, y*400) as box
        from (
            select
                st_rotate(
                    st_chaikinsmoothing(
                        st_makeline(array[
                            st_point(0,200),
                            st_point(100,300),
                            st_point(300,100),
                            st_point(400,200)]
                        ), 3, false),
                        pi()/2,
                        st_point(200, 200)
                ) as geom,
                st_envelope(st_collect(array[st_point(0,0), st_point(400,400)])) as box
        ) s,
        generate_series(0,4,1) y
    ),
    row as (
        /* generate a row*/
        select
            x, y,
            st_translate(geom, x*400, 0) as geom,
            st_translate(box, x*400, 0) as box
        from
            col,
            generate_series(0,15,1) x
    )
select
    /* ...and calculate data*/
    x, y, v, i, geom, st_numpoints(geom),
    c.avg_segment_len, lbl, box, p.box_anchor
from (
    select
        row.x, row.y, tolerance.v, tolerance.i,
        case
            when y=0 then st_reduceprecision(row.geom, tolerance.v)
            when y=1 then st_simplifypreservetopology(row.geom, tolerance.v)
            when y=2 then st_simplify(row.geom, tolerance.v)
            when y=3 then st_simplifyvw(row.geom, tolerance.v)
            else row.geom
        end as geom,
        case
            when y=0 and x = 0 then 'st_reduceprecision'
            when y=1 and x = 0 then 'st_simplifypreservetopology'
            when y=2 and x = 0 then 'st_simplify'
            when y=3 and x = 0 then 'st_simplifyvw'
            when x=0 then 'original'
        end as lbl,
        box
    from
        row
            left join lateral
                unnest(
                    array[
                        0.001, 0.01, 0.1, 0.5,
                        1, 5,
                        10, 20, 30, 40, 50,
                        100, 200, 300, 400, 500
                    ]
                ) with ordinality tolerance(v, i) on true
    where
        row.x = tolerance.i-1
) x
    join lateral (
        select
            avg(st_length((segs).geom)) as avg_segment_len    
        from
            st_dumpsegments(x.geom) segs
    ) c on true
    join lateral (
        select
            st_point(
                max(st_x(pts.geom)),
                min(st_y(pts.geom))
            ) as box_anchor
        from
            st_dumppoints(x.box) pts
    ) p on true



;