wm

Wang–Müller line generalization algorithm in PostGIS
Log | Files | Refs | README | LICENSE

visuals.sql (3805B) - Raw


      1 \set ON_ERROR_STOP on
      2 SET plpgsql.extra_errors TO 'all';
      3 
      4 -- wm_bbox clips a geometry by a bounding box around a given object,
      5 -- matching dimensions of A-class paper (1 by sqrt(2)).
      6 drop function if exists wm_bbox;
      7 create function wm_bbox(
      8   center geometry,
      9   scaledwidth float
     10 ) returns geometry as $$
     11 declare
     12   halfX float;
     13   halfY float;
     14 begin
     15   halfX = scaledwidth / 2;
     16   halfY = halfX * sqrt(2);
     17   return st_envelope(
     18     st_union(
     19       st_translate(center, halfX, halfY),
     20       st_translate(center, -halfX, -halfY)
     21     )
     22   );
     23 end $$ language plpgsql;
     24 
     25 -- wm_quadrant divides the given geometry to 4 rectangles
     26 -- and returns the requested quadrant following cartesian
     27 -- convention:
     28 --  +----------+
     29 --  | II  | I  |
     30 --- +----------+
     31 --  | III | IV |
     32 --  +-----+----+
     33 -- matching dimensions of A-class paper (1 by sqrt(2).
     34 drop function if exists wm_quadrant;
     35 create function wm_quadrant(
     36   geom geometry,
     37   quadrant integer
     38 ) returns geometry as $$
     39 declare
     40   xmin float;
     41   xmax float;
     42   ymin float;
     43   ymax float;
     44 begin
     45   xmin = st_xmin(geom);
     46   xmax = st_xmax(geom);
     47   ymin = st_ymin(geom);
     48   ymax = st_ymax(geom);
     49 
     50   if quadrant = 1 or quadrant = 2 then
     51     ymin = (ymin + ymax)/2;
     52   else
     53     ymax = (ymin + ymax)/2;
     54   end if;
     55 
     56   if quadrant = 2 or quadrant = 3 then
     57     xmax = (xmin + xmax)/2;
     58   else
     59     xmin = (xmin + xmax)/2;
     60   end if;
     61 
     62   return st_intersection(
     63     geom,
     64     st_makeenvelope(xmin, ymin, xmax, ymax, st_srid(geom))
     65   );
     66 end $$ language plpgsql;
     67 
     68 
     69 drop function if exists wm_salvisbbox;
     70 create function wm_salvisbbox(
     71   geom geometry,
     72   scaledwidth float
     73 ) returns geometry as $$
     74 declare
     75   ret geometry;
     76 begin
     77   with multismall as (
     78     select st_intersection(
     79       st_union(geom),
     80       wm_bbox(
     81         st_closestpoint(
     82           (select way from wm_rivers where name='Šalčia'),
     83           (select way from wm_rivers where name='Visinčia')
     84         ),
     85         scaledwidth
     86       )
     87     ) ways
     88   )
     89   -- protecting against very small bends that were cut
     90   -- in the corner of the picture
     91   select st_union(a.geom)
     92   from st_dump((select ways from multismall)) a
     93   where st_length(a.geom) >= 100
     94   into ret;
     95 
     96   return ret;
     97 end $$ language plpgsql;
     98 
     99 delete from wm_debug where name like 'salvis%';
    100 delete from wm_visuals where name like 'salvis%';
    101 insert into wm_visuals(name, way) values
    102   ('salvis-grpk10', (
    103     wm_salvisbbox(
    104       (select st_union(way) from wm_rivers where name in ('Šalčia', 'Visinčia')),
    105       :scaledwidth
    106     )
    107   )),
    108   ('salvis-grpk50', (
    109     wm_salvisbbox(
    110       (select st_union(way) from wm_rivers_50 where name in ('Šalčia', 'Visinčia')),
    111       :scaledwidth
    112     )
    113   )),
    114   ('salvis-grpk250', (
    115     wm_salvisbbox(
    116       (select st_union(way) from wm_rivers_250 where name in ('Šalčia', 'Visinčia')),
    117       :scaledwidth
    118     )
    119   ));
    120 
    121 do $$
    122 declare
    123   i integer;
    124   geom1 geometry;
    125   geom2 geometry;
    126   geom3 geometry;
    127 begin
    128   foreach i in array array[16, 32, 64, 256] loop
    129     geom1 = st_simplify((select way from wm_visuals where name='salvis-grpk10'), i);
    130     geom2 = st_simplifyvw((select way from wm_visuals where name='salvis-grpk10'), i*i);
    131     insert into wm_visuals(name, way) values
    132       ('salvis-dp'        || i, geom1),
    133       ('salvis-dpchaikin' || i, st_chaikinsmoothing(geom1, 5)),
    134       ('salvis-vw'        || i, geom2),
    135       ('salvis-vwchaikin' || i, st_chaikinsmoothing(geom2, 5));
    136   end loop;
    137 
    138   -- more than 220 doesn't work, because there is an exaggerated bend near
    139   -- Šalčia-Visinčia crossing, and it "exaggerates" to the
    140   -- other river.
    141   foreach i in array array[75, 220] loop
    142     geom3 = st_simplifywm((select way from wm_visuals where name='salvis-grpk10'), i, null, 'salvis-wm' || i);
    143     insert into wm_visuals(name, way) values
    144       ('salvis-wm'          || i, geom3);
    145   end loop;
    146 end $$ language plpgsql;