-- Add some input variables
-- Epsilon is the merge distance of two close points
set @EPSILON = 1;
-- River width is the minimal river width to take account for
set @MINIMAL_RIVER_WIDTH = 4;
-- Centrality is a network importance coefficient filtering between 0 and 1.
set @RATIO_CENTRALITY = 0.8;
-- Merge adjacent polygons, the simplify them and final densify their contours. Densification is done in order to create triangles where sides must fit into the width 
of the polygons. A theoretical density of 0m will create a valid skeleton of polygons.
drop table if exists watersurf_accum;
create table watersurf_accum as select ST_SIMPLIFYPRESERVETOPOLOGY(ST_BUFFER(ST_UNION(ST_ACCUM(the_geom)), 0),@EPSILON) the_geom from POLYGONS;
-- Create voronoi edges
drop table if exists water_voronoi;
create table water_voronoi as select 
ST_VORONOI(ST_DELAUNAY(ST_UPDATEZ(ST_REMOVEREPEATEDPOINTS(ST_PRECISIONREDUCER(ST_TOMULTIPOINT(ST_DENSIFY(ST_ACCUM(THE_GEOM),@MINIMAL_RIVER_WIDTH / 2)),@EPSILON)),0)), 
1) the_geom from watersurf;
-- Explode the MultiLineStrings into multiple rows of LineString
drop table if exists water_voronoi_edges;
create table water_voronoi_edges as select * from st_explode('water_voronoi');
-- Create optimisation structure in order to speedup filtering
SET @EXTENT = SELECT ST_EXTENT(THE_GEOM) THE_GEOM FROM watersurf_accum;
SET @WIDTH = SELECT ST_XMAX(@EXTENT) - ST_XMIN(@EXTENT);
SET @HEIGHT = SELECT ST_YMAX(@EXTENT) - ST_YMIN(@EXTENT);
SET @CELL_SIZE = GREATEST(@WIDTH / 30, @HEIGHT/30);
drop table if exists GRIDED_POLYGON;
CREATE TABLE GRIDED_POLYGON AS SELECT ST_INTERSECTION(G.THE_GEOM, P.THE_GEOM) THE_GEOM FROM WATERSURF_ACCUM P, ST_MAKEGRID(@EXTENT,@CELL_SIZE, @CELL_SIZE) G WHERE 
ST_INTERSECTS(G.THE_GEOM, P.THE_GEOM);
create spatial index on GRIDED_POLYGON(THE_GEOM);
-- Keep only the voronoi edges that are inside the polygons
drop table if exists water_skel;
create table water_skel(THE_GEOM GEOMETRY) as select V.THE_GEOM THE_GEOM FROM water_voronoi_edges V WHERE ST_CONTAINS((SELECT ST_UNION(ST_ACCUM(W.THE_GEOM)) the_geom 
FROM GRIDED_POLYGON W WHERE v.the_geom && w.the_geom), v.the_geom);
-- Union and simplify skeleton
drop table if exists water_skel_simple_accum;
create table water_skel_simple_accum as select ST_LINEMERGE(ST_ACCUM(THE_GEOM)) THE_GEOM FROM water_skel;
drop table if exists water_skel_simple;
create table water_skel_simple(PK SERIAL PRIMARY KEY, THE_GEOM GEOMETRY) as select null, the_geom from ST_EXPLODE('water_skel_simple_accum');
drop table if exists water_skel_simple_accum;
drop table if exists graph_waternetwork;
drop table if exists WATER_SKEL_SIMPLE_EDGES, WATER_SKEL_SIMPLE_NODES;
call ST_GRAPH('WATER_SKEL_SIMPLE');
drop table if exists EDGES_WEIGHT;
create table edges_weight as select E.*, ST_LENGTH(THE_GEOM) weight FROM WATER_SKEL_SIMPLE_EDGES E, WATER_SKEL_SIMPLE W WHERE W.PK = E.EDGE_ID ;
-- This call is time and memory consuming
drop table if exists EDGES_WEIGHT_EDGE_CENT,EDGES_WEIGHT_NODE_CENT;
call ST_GraphAnalysis('EDGES_WEIGHT', 'undirected', 'weight');
drop table if exists central_skel;
create table central_skel as select the_geom, LEAST(ENDNODE.BETWEENNESS, STARTNODE.BETWEENNESS) BETWEENNESS from EDGES_WEIGHT_NODE_CENT ENDNODE,EDGES_WEIGHT_NODE_CENT 
STARTNODE,WATER_SKEL_SIMPLE_EDGES E, WATER_SKEL_SIMPLE W WHERE E.START_NODE = STARTNODE.NODE_ID AND E.END_NODE = ENDNODE.NODE_ID AND W.PK = E.EDGE_ID AND 
LEAST(ENDNODE.BETWEENNESS, STARTNODE.BETWEENNESS) > @RATIO_CENTRALITY;