subdividing "empty" space by expanding polygons in PostGIS

expanding-polygons.sql

/** Subdividing "emtpy" space by expanding polygons
*   -----------------------------------------------
*   This is the code-companion that is used in 
*   https://tkardi.ee/writeup/post/2018/07/21/subdividing-space/
*
*   Input data includes shapefiles from the Estonian Land Board's homepage 
*   (downloaded May 2018): 
*   - Estonian third level administrative units (settlements) with coastline 
*     borders: "asustusüksus", available for download from 
*     https://geoportaal.maaamet.ee/docs/haldus_asustus/asustusyksus_shp.zip
*   - Estonian second level administrative units (municipalities) with outer
*     boundaries expanded to the sea: "omavalitsus", available for download 
*     from
*     https://geoportaal.maaamet.ee/docs/haldus_asustus/omavalitsus_valispiirini_shp.zip
*   The datasets do not seem to have a clear indication of a license, but the 
*   download page states the data is free to use as long as it's attributed 
*   to the Land Board and data checkout date is given. More information on 
*   these datasets is available at 
*   https://geoportaal.maaamet.ee/est/Andmed-ja-kaardid/Haldus-ja-asustusjaotus-p119.html
*   (in Estonian). Unfortunately the English version of the page seems to be 
*   a bit out of date: 
*   https://geoportaal.maaamet.ee/eng/Maps-and-Data/Administrative-and-Settlement-Division-p312.html
*
*   The following SQL code is published under the Unlicense.
*/

/* SPDX-License-Identifier: Unlicense */


/* extract rings of settlement areas */

drop table if exists ehak.ay_rings;
create table ehak.ay_rings as
select
    row_number() over()::int as oid, bar.*
from (
    select
        okood, mkood, array_agg(akood) as akoods,
        (st_dump(st_node(st_collect(geom)))).geom
    from (
        select
            akood, okood, mkood,
            st_exteriorring((st_dumprings((st_dump(geom)).geom)).geom) as geom
        from ehak.asustusyksus
    ) foo
    group by okood, mkood
) bar;

-- create index
create index sidx__ay_rings on ehak.ay_rings using gist(geom);


/* create a "mesh" of settlements (i.e. a single line border
between two neighbouring polygons) */

drop table if exists ehak.ay_rings_mesh;
create table ehak.ay_rings_mesh as
select row_number() over()::int as oid, bar.*
from (
    select
        (st_dump(st_linemerge(st_collect(geom)))).geom as geom
    from (
        select
            (st_dump(st_union(geom))).geom as geom
        from ehak.ay_rings
    ) foo
) bar;

-- create index
create index sidx__ay_rings_mesh on ehak.ay_rings_mesh using gist(geom);


/* create half-way points, and find their territorial belonging */

drop table if exists ehak.ay_rings_mesh_points;
create table ehak.ay_rings_mesh_points as
select oid, st_lineinterpolatepoint(geom, 0.5) as geom
from ehak.ay_rings_mesh;

create index sidx__ay_rings_mesh_points on
    ehak.ay_rings_mesh_points using gist(geom);

alter table ehak.ay_rings_mesh_points add column akoods varchar[];
alter table ehak.ay_rings_mesh_points add column okoods varchar[];
alter table ehak.ay_rings_mesh_points add column count integer;

alter table ehak.ay_rings_mesh_points
    add constraint pk__ay_rings_mesh_points primary key (oid);

update ehak.ay_rings_mesh_points set
    okoods = n.okoods,
    akoods = n.akoods,
    count = n.count
from (
    select
        p.oid, array_agg(n.okood) as okoods, array_agg(n.akood) as akoods,
        count(1)
    from ehak.asustusyksus n, ehak.ay_rings_mesh_points p
    where st_intersects(st_buffer(p.geom, 0.0001), n.geom)
    group by p.oid
) n
where n.oid = ay_rings_mesh_points.oid;


/* find vertices we're going to do Voronoi on */

drop table if exists ehak.ay_expand_points;
create table ehak.ay_expand_points as
select
    akood, okood,
    (st_dumppoints(
        st_linesubstring(geom, 0.5/st_length(geom), 1-0.5/st_length(geom))
    )).*
from (
    select
        m.geom, p.akoods[1] as akood, p.okoods[1] as okood
    from ehak.ay_rings_mesh_points p, ehak.ay_rings_mesh m
    where p.oid = m.oid and p.count = 1
) outs;

-- add index on geometry
create index sidx__ay_expand_points on ehak.ay_expand_points using gist (geom);


/* create Voronoi polygons for those and trim them to shape */

drop table if exists ehak.ay_expand_voros;
create table ehak.ay_expand_voros as
select
    row_number() over()::int as oid, baz.okood, pts.akood,
    (st_dump(
        st_intersection(baz.geom, baz.cutter)
    )).geom
from (
    select bar.okood,
        (st_dump(
            st_voronoipolygons(bar.geom, 0.0, q.geom)
        )).geom, q.geom as cutter
    from (
        select okood, st_collect(st_snaptogrid(foo.geom, 0.001)) as geom
        from ehak.ay_expand_points foo
        group by okood
    ) bar, ehak.omavalitsus q
    where bar.okood = q.okood
) baz, ehak.ay_expand_points pts
where
    pts.okood = baz.okood and
    st_within(pts.geom, baz.geom);

-- add geometry index
create index sidx__ay_expand_voros on ehak.ay_expand_voros using gist (geom);


/* union them together based on the settlement identifier */

drop table if exists ehak.ay_expand_unions;
create table ehak.ay_expand_unions as
select akood, okood, st_union(geom) as geom
from ehak.ay_expand_voros
group by akood, okood;


/* preparations for cookie-cutting */

alter table ehak.ay_expand_unions add column cutter geometry(Geometry, 3301);
alter table ehak.ay_expand_unions add column self geometry(Geometry, 3301);
alter table ehak.ay_expand_unions add column unioned geometry(Geometry, 3301);

-- self, which is the original geometry of the settlement
update ehak.ay_expand_unions set
    self = n.geom
from (
    select akood, st_union(geom) as geom
    from ehak.asustusyksus
    group by akood
) n
where ay_expand_unions.akood = n.akood;


-- the cutter, which are all other settlement areas in this municipality
update ehak.ay_expand_unions set
    cutter = foo.geom
from (
    select s.akood, st_union(o.geom) as geom
    from
        ehak.asustusyksus s,
        ehak.asustusyksus o,
        (select distinct akood from ehak.ay_expand_unions) n
    where
        s.akood != o.akood and
        s.okood = o.okood and
        s.akood = n.akood
    group by s.akood
) foo
where foo.akood = ay_expand_unions.akood;

-- and unioned, which is the point set union of the original settlement polygon
-- and it's expanded counterpart  

update ehak.ay_expand_unions set
    unioned = n.geom
from (
    select akood, st_union(geom) as geom
    from (
        select akood, (st_dump(st_union(geom, self))).geom
        from ehak.ay_expand_unions
    ) foo
    where geometrytype(geom) = 'POLYGON'
    group by akood
) n
where ay_expand_unions.akood = n.akood;


/* and finally cut everything to shape (cookie-cutting) */

drop table if exists ehak.ay_expand_trimmed;
create table ehak.ay_expand_trimmed as
select akood, okood, st_union(geom) as geom
from (
    select akood, okood,
    case
        when geom is not null and cutter is not null then
            (st_dump(st_difference(geom, cutter))).geom
        when geom is not null and cutter is null then
            geom
        else null
    end as geom
    from (
        select akood, okood, unioned as geom,
        cutter as cutter
        from ehak.ay_expand_unions) a
    ) foo
where geometrytype(geom) = 'POLYGON' or geometrytype(geom) = 'MULTIPOLYGON'
group by akood, okood;

unlicense

This is free and unencumbered software released into the public domain.

Anyone is free to copy, modify, publish, use, compile, sell, or
distribute this software, either in source code form or as a compiled
binary, for any purpose, commercial or non-commercial, and by any
means.

In jurisdictions that recognize copyright laws, the author or authors
of this software dedicate any and all copyright interest in the
software to the public domain. We make this dedication for the benefit
of the public at large and to the detriment of our heirs and
successors. We intend this dedication to be an overt act of
relinquishment in perpetuity of all present and future rights to this
software under copyright law.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
OTHER DEALINGS IN THE SOFTWARE.

For more information, please refer to <http://unlicense.org/>