Skip to content

Instantly share code, notes, and snippets.

@tkardi
Last active July 30, 2018 10:27
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save tkardi/4d09093ec823ee44c5f35c56ae223406 to your computer and use it in GitHub Desktop.
Save tkardi/4d09093ec823ee44c5f35c56ae223406 to your computer and use it in GitHub Desktop.
Calculate a polygon "mesh" with left/right side properties in PostGIS
/** Calculate a polygon "mesh" with left/right side properties in PostGIS
* ----------------------------------------------------------------------
* Based on Estonian third level administrative units data downloaded
* from the Estonian Land Board's homepage (May 2018) @
* https://geoportaal.maaamet.ee/eng/Maps-and-Data/Administrative-and-Settlement-Division-p312.html
* Direct link to the shapefile - https://geoportaal.maaamet.ee/docs/haldus_asustus/asustusyksus_shp.zip
*
* The following SQL code is published under the Unlicense.
*/
/* SPDX-License-Identifier: Unlicense */
-- polygons to linestrings
drop table if exists ay_aslines;
create table ay_aslines as
select st_exteriorring((st_dumprings(shape)).geom) as geom
from asustusyksus;
-- merge all lines, dump to singleparts
-- yields 2-node linestrings
drop table if exists ay_lines_merged;
create table ay_lines_merged as
select (st_dump(st_union(geom))).geom
from ay_aslines;
-- sow 2-node linestrings together
-- yield noded singlepart linestrings
drop table if exists ay_lines_merged_full;
create table ay_lines_merged_full as
select
path[1] as id, geom
from (
select (st_dump(st_linemerge(st_collect(geom)))).*
from ay_lines_merged
) foo;
-- calculate sidedness
-- shift left
drop table if exists ay_lines_shift_left;
create table ay_lines_shift_left as
select id, (st_dump(st_offsetcurve(geom, 0.01))).geom as geom
from ay_lines_merged_full;
-- ... and shift right
drop table if exists ay_lines_shift_right;
create table ay_lines_shift_right as
select id, (st_dump(st_offsetcurve(geom, -0.01))).geom as geom
from ay_lines_merged_full;
-- columns for holding the midpoint of shifted lines
alter table ay_lines_shift_left add column midpoint geometry(Point, 3301);
alter table ay_lines_shift_right add column midpoint geometry(Point, 3301);
-- calculate midpoints on the left
update ay_lines_shift_left set
midpoint = st_lineinterpolatepoint(geom, 0.5);
-- .. and on the right
update ay_lines_shift_right set
midpoint = st_lineinterpolatepoint(geom, 0.5);
-- columns to hold the sidedness info
alter table ay_lines_shift_left add column a3_id varchar(4);
alter table ay_lines_shift_right add column a3_id varchar(4);
create index sidx__ay_lines_shift_left__midpoint
on ay_lines_shift_left using gist (midpoint);
create index sidx__ay_lines_shift_right__midpoint
on ay_lines_shift_right using gist (midpoint);
-- find the polygons that the shifted border midpoints are within
-- on the left
update ay_lines_shift_left set
a3_id = ay.akood
from asustusyksus ay
where st_within(ay_lines_shift_left.midpoint, ay.shape);
-- .. and on the right
update ay_lines_shift_right set
a3_id = ay.akood
from asustusyksus ay
where st_within(ay_lines_shift_right.midpoint, ay.shape);
-- add some more columns to the output table
alter table ay_lines_merged_full add column left_a3 varchar(4);
alter table ay_lines_merged_full add column right_a3 varchar(4);
alter table ay_lines_merged_full add column left_a2 varchar(4);
alter table ay_lines_merged_full add column right_a2 varchar(4);
alter table ay_lines_merged_full add column left_a1 varchar(4);
alter table ay_lines_merged_full add column right_a1 varchar(4);
-- ... and update their values for A3, A2, A1 admin levels
update ay_lines_merged_full set
left_a3 = n.a3_id
from ay_lines_shift_left n
where n.id = ay_lines_merged_full.id;
update ay_lines_merged_full set
right_a3 = n.a3_id
from ay_lines_shift_right n
where n.id = ay_lines_merged_full.id;
update ay_lines_merged_full set
left_a2 = kov.okood
from (select distinct akood, okood from asustusyksus) kov
where kov.akood = ay_lines_merged_full.left_a3
update ay_lines_merged_full set
right_a2 = kov.okood
from (select distinct akood, okood from asustusyksus) kov
where kov.akood = ay_lines_merged_full.right_a3
update ay_lines_merged_full set
left_a1 = mk.mkood
from (select distinct akood, mkood from ay) mk
where mk.akood = ay_lines_merged_full.left_a3;
update ay_lines_merged_full set
right_a1 = mk.mkood
from (select distinct akood, mkood from ay) mk
where mk.akood = ay_lines_merged_full.right_a3;
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/>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment