back to start page [[..]] [[PageOutline]] * [#PROBLEM Fixing a problem with the database] * [#FIRST First trial with two layers] * [#PFAF Trials for the pfafsteter system] * [#ROUTING pg_routing] * [#LTREE ltree] * [#BASSIN Select some bassins to calculate topology] * [#PORTUGAL Trial with Portugal layers] === Tickets === #90 #97 = Import layer postgis = == first fixing a problem with the database == [=#PROBLEM] At first the following failed and I had to fix a problem in the database ~~ spatial_ref_sys is empty in the eda.0 database => collecting it from another database ~~ {{{#!sh pg_dump -U postgres -f "spatial_ref_sys.sql" --data-only --table public.spatial_ref_sys --inserts --verbose eda2.2 psql -U postgres -f spatial_ref_sys.sql eda2.0 }}} this does not work, it's a known bug https://trac.osgeo.org/postgis/ticket/1815 I'm using the sql file from the postgis installation. For this we will try on a subset in the oria {{{#!sh psql -U postgres -f "C:/Program Files/PostgreSQL/9.6/share/contrib/postgis-2.4/spatial_ref_sys.sql" eda2.0 }}} == First trial with two layers == [=#FIRST] Our objective is to test the best way to route the river network {{{#!sql -- first create a new schema for spain psql -U postgres -c "create schema spain" eda2.0 -- shp2pgsql with create index option -- Subcuencas de ríos completos clasificadas según Pfafstetter modificado shp2pgsql -s 4258 -I A_cuencas_rios_Atl_Norte.shp spain.a_cuencas_rios_atl_norte > a_cuencas_rios_atl_norte.sql psql -U postgres -f a_cuencas_rios_atl_norte.sql eda2.0 -- Tramos de ríos de España clasificados según Pfafstetter modificado shp2pgsql -s 4258 -W LATIN1 -I A_Rios_v2.shp spain.a_rios_v2 > a_rios_v2.sql psql -U postgres -f a_rios_v2.sql eda2.0 }}} in fact the system of projection indicated on the inspire layers is wrong it is 25830 http://www.mapama.gob.es/es/cartografia-y-sig/ide/directorio_datos_servicios/caracteristicas_wms.aspx {{{#!sql select UpdateGeometrySRID('spain','a_cuencas_rios_atl_norte', 'geom', 25830 ); select UpdateGeometrySRID('spain','a_rios_v2', 'geom', 25830 ); }}} Now we need to select only the Oria and neighbouring rivers for our trial {{{#!sql set search_path to public, spain ;-- using this I will not have to say that the tables are in the spain scheme select * from a_cuencas_rios_atl_norte limit 10 select * from a_cuencas_rios_atl_norte where nom_rio_1 like '%ORIA%'; }}} ||802F802||"A"||172988||"1003714"||"ORIA IBAIA"||""||2707||2707||"A"||175359||"10037294"||...|| 10037 corresponds to the Biscay coast of Spain, 1003714 correspond à Oria {{{#!sql select * from a_cuencas_rios_atl_norte where pfafrio like '10037%'; select * from a_cuencas_rios_atl_norte where pfafrio like '1003714%'; }}} Creating a table corresponding to ORIA Table of basins {{{#!sql drop table if exists spain.oria_cuenca; create table spain.oria_cuenca as select * from a_cuencas_rios_atl_norte where pfafrio like '1003714%';--422 CREATE INDEX oria_cuenca_geom_idx ON oria_cuenca USING gist (geom); }}} Table of rivers {{{#!sql drop table if exists spain.oria_a_rios_v2; create table spain.oria_a_rios_v2 as select * from a_rios_v2 where pfafrio like '1003714%';--422 CREATE INDEX oria_a_rios_v2_geom_idx ON oria_a_rios_v2 USING gist (geom); }}} [[Image(source:eda/data/Docs/trac/sudoang/oria.png,200px)]] Creating a table corresponding to UROLA Table of basins {{{#!sql create table spain.urola_cuenca as select * from a_cuencas_rios_atl_norte where pfafrio like '1003716%'; CREATE INDEX urola_cuenca_geom_idx ON urola_cuenca USING gist (geom); }}} Table of rivers {{{#!sql create table spain.urola_a_rios_v2 as select * from a_rios_v2 where pfafrio like '1003716%'; CREATE INDEX urola_a_rios_v2_geom_idx ON urola_a_rios_v2 USING gist (geom); }}} [[Image(source:eda/data/Docs/trac/sudoang/urola.png,200px)]] == Trials for the pfafsteter system == [=#PFAF] To be able to explore the topology we could use pfafstetter {{{#!sql /*----------------------------------------- Function to get all segments downstream from one one in the same subbasin Problem when the subbasin changes it is difficult to go downstream -----------------------------------------*/ CREATE OR REPLACE FUNCTION downstream_segments(_id text) RETURNS TABLE(pfafcuen character varying(254), lngtramo_m numeric) AS $$ BEGIN RAISE NOTICE 'downstream segment for %', _id; -- prints the selected segment RETURN QUERY -- RETURN QUERY allows to return a table without loop SELECT r.pfafcuen, r.lngtramo_m from oria_a_rios_v2 r WHERE substring(_id::text,1,char_length(_id::text)-1)= substring(r.pfafcuen::text,1,char_length(r.pfafcuen::text)-1) -- the selection is done on a subset eg. 10011 10013 10015 10017 will be selected using string 1001 AND _id>r.pfafcuen -- we select only segments downstream from the current one AND (r.pfafcuen::numeric % 2) = 1 -- even number order by r.pfafcuen::numeric DESC; END $$ LANGUAGE plpgsql VOLATILE; select * from downstream_segments('10037146457'); }}} this query returns the three odd nodes dowstream from '10037146457' ||pfafcuen||lngtramo_m|| ||10037146451||747|| ||10037146453||304|| ||10037146455||437|| Problem from that step we need to go to another number 10037146451 => 1003714643 {{{#!sql -- this query returns nothing with d1 as (SELECT pfafcuen, lngtramo_m, (substring(pfafcuen::text,1,char_length(pfafcuen::text)-1)::integer-2)::character varying(254) as pfafcuen_downstream from downstream_segments('10037146457')) select * from downstream_segments((select min(pfafcuen_downstream) FROM d1 limit 1)); --because this query returns nothing select downstream_segments('1003714643'); }}} Before trying too much the following problem arises [[BR]] where `1037146(47)` is immediately upstream from `1037146(457)` [[BR]] [[Image(source:eda/data/Docs/trac/sudoang/network_topology.PNG,400px)]][[BR]] == pg_routing == [=#ROUTING] https://docs.pgrouting.org/2.5/en/pgRouting-concepts.html#description-of-the-edges-sql-query-for-dijkstra-like-functions {{{#!sql ------------------------------------- -- pgrouting installation and configuration on table CREATE EXTENSION pgrouting; ALTER TABLE spain.oria_a_rios_v2 DROP COLUMN source; ALTER TABLE spain.oria_a_rios_v2 DROP COLUMN target; ALTER TABLE spain.oria_a_rios_v2 ADD COLUMN source integer; ALTER TABLE spain.oria_a_rios_v2 ADD COLUMN target integer; CREATE INDEX oria_a_rios_v2_source_idx ON spain.oria_a_rios_v2 (source); CREATE INDEX oria_a_rios_v2_target_idx ON spain.oria_a_rios_v2 (target); --------------------------------------------------------------------- -- createTopology also creates spain.oria_a_rios_v2_vertices_pgr ------------------------------------------------------------------ set search_path to spain,public; SELECT pgr_createTopology('oria_a_rios_v2', 0.0001, 'geom', 'gid');--1.4s --------------------------------------------------------------------- -- routing this is done using nodes, FALSE is necessery otherwise fails (using a directed path= ------------------------------------------------------------------ SELECT seq, node, edge, agg_cost, gid, pfafcuen FROM pgr_dijkstra( 'SELECT gid as id, source, target, st_length(geom) as cost FROM oria_a_rios_v2', 5, 540, FALSE ) pt JOIN oria_a_rios_v2 r ON pt.edge = r.gid order by seq; }}} ||"seq"||"node"||"edge"||"agg_cost"||"gid"||"pfafcuen"|| ||1||5||774||0||774||"1003714113"|| ||2||28||783||340.239109235141||783||"1003714115"|| ||3||33||916||434.71909815275||916||"1003714117"|| ||4||61||920||1381.42299559645||920||"100371413"|| ||5||62||956||1931.41349601216||956||"10037141511"|| ||6||54||992||2751.07786782504||992||"10037141513"|| ||7||53||1023||3006.69704002442||1023||"10037141515"|| ||8||16||1028||3923.51604452287||1028||"10037141517"|| ||9||15||962||4550.98432547035||962||"10037141531"|| ||10||47||879||5241.61915518433||879||"10037141533"|| ||11||39||880||5777.24646955956||880||"10037141535"|| == ltree == [=#LTREE] {{{#!sql CREATE EXTENSION ltree; ALTER TABLE spain.oria_a_rios_v2 ADD COLUMN chemin ltree; CREATE INDEX chemin_gist_oria_idx ON spain.oria_a_rios_v2 USING GIST(chemin); CREATE INDEX chemin_oria_idx ON spain.oria_a_rios_v2 USING btree(chemin); -- WHAT IS THE MOST DOWNSTREAM SEGMENT ? -- reprojection of the wise layer -- still a bit far create view european_wise2008.coastal_waters_25830 as select gid, code, name, st_transform(the_geom,25830) as geom from european_wise2008.costal_waters where cty_id='ES'; st_distance(c.geom, r.geom) -- renvoit l'id du point le plus proche de la mer select id from ( select v.id, st_distance(c.geom, v.the_geom) as dist from european_wise2008.coastal_waters_25830 c JOIN spain.oria_a_rios_v2_vertices_pgr v ON ST_DWithin(c.geom, v.the_geom, 1000)) jointure order by dist limit 1 --calcule le nombre de lignes de la table Select count(*) from spain.oria_a_rios_v2 ; select * from spain.oria_a_rios_v2 limit 10; select * from spain.oria_a_rios_v2 where target=6 -- this query returns the ltree "1003714113.1003714115.1003714117(...)" SELECT text2ltree(string_agg(pfafcuen, '.')) AS pfafcuenlist FROM (SELECT pfafcuen FROM pgr_dijkstra( 'SELECT gid as id, source, target, st_length(geom) as cost FROM spain.oria_a_rios_v2', 5, 540, FALSE ) pt JOIN spain.oria_a_rios_v2 r ON pt.edge = r.gid order by seq) sub -- function to create ltree from two nodes from and to DROP function get_path(integer,integer); CREATE OR REPLACE FUNCTION spain.get_path(_from integer, _to integer) RETURNS SETOF ltree AS $$ BEGIN RETURN QUERY SELECT text2ltree(string_agg(pfafcuen, '.')) AS pfafcuenlist FROM (SELECT pfafcuen FROM pgr_dijkstra( 'SELECT gid as id, source, target, st_length(geom) as cost FROM spain.oria_a_rios_v2', _from, _to, FALSE ) pt JOIN spain.oria_a_rios_v2 r ON pt.edge = r.gid order by seq) sub; END $$ LANGUAGE plpgsql VOLATILE; select spain.get_path(5,540); SELECT seq, node, edge, agg_cost, gid, pfafcuen FROM pgr_dijkstra( 'SELECT gid as id, source, target, st_length(geom) as cost FROM spain.oria_a_rios_v2', 5, 540, FALSE ) pt JOIN spain.oria_a_rios_v2 r ON pt.edge = r.gid order by seq; UPDATE spain.oria_a_rios_v2 set chemin=get_path(6,540) where source=540 /* didn't work but keep it for later checking -- rio the most downstream with pfaf_length as ( select length(pfafrio) l,pfafrio from spain.oria_a_rios_v2) select pfafrio from pfaf_length order by l limit 1 -- removing the rio with shortened_pfaf as ( select regexp_replace(pfafcuen,'1003714','') pfaf_short from spain.oria_a_rios_v2) select substring(pfaf_short,1,3) from shortened_pfaf select substring (1,2, */ select format('Select count(*) FROM %1$I','toto') /* TODO test that the routing algorythm always set the most downstream point as a target ... If so develop a formal testing within the function */ select spain.get_path(6,source) from spain.oria_a_rios_v2 limit 10 }}} {{{ "get_path" "1003714111.1003714113.10037141141.10037141142" "1003714111.1003714113.1003714115.1003714117.100371413.10037141411.10037141413.10037141421.10037141422" "1003714111.1003714113.1003714115.1003714117.100371413.10037141511.10037141513.10037141515.10037141517.10037141531.10037141533.10037141535.10037141537.10037141539.10037141551.10037141553.10037141555.10037141557.10037141571.10037141573.10037141575.1003714157 (...)" "1003714111.1003714113.1003714115.1003714117.100371413.10037141511.10037141513.10037141515.10037141517.10037141531.10037141533.10037141535.10037141537.10037141539.10037141551.10037141553.10037141555.10037141557.10037141571.10037141573.10037141575.1003714157 (...)" "1003714111.1003714113.1003714115.1003714117.100371413.10037141511.10037141513.10037141515.10037141517.10037141531.10037141533.10037141535.10037141537.10037141539.10037141551.10037141553.10037141555.10037141557.10037141571.10037141573.10037141575.1003714157 (...)" "1003714111.1003714113.1003714115.1003714117.100371413.10037141511.10037141513.10037141515.10037141517.10037141531.10037141533.10037141535.10037141537.10037141539.10037141551.10037141553.10037141555.10037141557.10037141571.10037141573.10037141575.1003714157 (...)" "1003714111.1003714113.1003714115.1003714117.100371413.10037141511.10037141513.10037141515.10037141517.10037141531.10037141533.10037141535.10037141537.10037141539.10037141551.10037141553.10037141555.10037141557.10037141571.10037141573.10037141575.1003714157 (...)" "1003714111.1003714113.1003714115.1003714117.100371413.10037141511.10037141513.10037141515.10037141517.10037141531.10037141533.10037141535.10037141537.10037141539.10037141551.10037141553.10037141555.10037141557.10037141571.10037141573.10037141575.1003714157 (...)" "1003714111.1003714113.1003714115.1003714117.100371413.10037141511.10037141513.10037141515.10037141517.10037141531.10037141533.10037141535.10037141537.10037141539.10037141551.10037141553.10037141555.10037141557.10037141571.10037141573.10037141575.1003714157 (...)" "1003714111.1003714113.1003714115.1003714117.100371413.10037141511.10037141513.10037141515.10037141517.10037141531.10037141533.10037141535.10037141537.10037141539.10037141551.10037141553.10037141555.10037141557.10037141571.10037141573.10037141575.1003714157 (...)" }}} {{{ Query returned successfully: 805 rows affected, 03:17 minutes execution time. }}} == Select some bassins to calculate topology == [=#BASSIN] {{{#!sql set search_path to public, spain, european_wise2008 ; -- location of the tables in the scheme -- First trial with two layers ("a_cuencas_rios_atl_norte" and "a_rios_v2") -- Function to select some bassins to calculate topology -- Rivers --select * from a_rios_v2 limit 10; -- visualize 10 first rows --ALTER TABLE a_rios_v2 DROP COLUMN pfafbas; -- delete columns --ALTER TABLE a_rios_v2 DROP COLUMN rio; ALTER TABLE a_rios_v2 ADD COLUMN pfafbas varchar; ALTER TABLE a_rios_v2 ADD COLUMN rio varchar; UPDATE a_rios_v2 SET pfafbas = SUBSTR(pfafrio, 1, 5), -- New columns are created as character! rio = SUBSTR(pfafrio, 6, 2) --SELECT *, pfafrio, (SUBSTR(pfafrio, 1, 5) = '10037%') FROM a_rios_v2; -- create a column boolean to say TRUE/FALSE --SELECT *, SUBSTR(pfafrio, 1, 5) "pfafbas", SUBSTR(pfafrio, 6, 2) "rios" --FROM a_rios_v2; -- Test with '10037%' that corresponds to the Biscay coast of Spain: -- Table of Basins drop table if exists biscay_coast_cuenca; create table biscay_coast_cuenca as select * from a_cuencas_rios_atl_norte where pfafrio like '10037%'; CREATE INDEX biscay_coast_cuenca_geom_idx ON biscay_coast_cuenca USING gist (geom); -- Table of rivers drop table if exists biscay_coast_a_rios_v2; create table biscay_coast_a_rios_v2 as select * from a_rios_v2 where pfafrio like '10037%'; CREATE INDEX biscay_coast_a_rios_v2_geom_idx ON biscay_coast_a_rios_v2 USING gist (geom); }}} Check in QGis the rivers created according to the 6th and 7th position of code in pfafrio [[BR]] [[Image(source:eda/data/Docs/trac/sudoang/test_biscay_coast.png,400px)]] [[BR]] pgRouting, configuration on table and ltree {{{#!sh ALTER TABLE biscay_coast_a_rios_v2 DROP COLUMN if exists source; ALTER TABLE biscay_coast_a_rios_v2 DROP COLUMN if exists target; ALTER TABLE biscay_coast_a_rios_v2 ADD COLUMN source integer; ALTER TABLE biscay_coast_a_rios_v2 ADD COLUMN target integer; CREATE INDEX biscay_coast_a_rios_v2_source_idx ON biscay_coast_a_rios_v2 (source); CREATE INDEX biscay_coast_a_rios_v2_target_idx ON biscay_coast_a_rios_v2 (target); --------------------------------------------------------------------- -- createTopology also creates biscay_coast_a_rios_v2_vertices_pgr ------------------------------------------------------------------ SELECT pgr_createTopology('biscay_coast_a_rios_v2', 0.0001, 'geom', 'gid'); }}} Get error creating topology == Trial with Portugal layers == [=#PORTUGAL] In Command prompt(cmd) we need to create a new schema and import the shapefiles to Postgres {{{#!sql -- Create a new schema for Portugal psql -U postgres -c "create schema portugal" eda2.0 -- shp2pgsql with create index option -- Bacias SNIRH (to import shapefiles you must be in the directory where the files are found) shp2pgsql -s 3763 -I AtAgua_Bacias_bacias_snirh_PC.shp portugal.AtAgua_Bacias_bacias_snirh_PC > AtAgua_Bacias_bacias_snirh_PC.sql psql -U postgres -f AtAgua_Bacias_bacias_snirh_PC.sql eda2.0 -- Rede hidrografica GeoCodificada shp2pgsql -s 3763 -W LATIN1 -I netElementL.shp portugal.netElementL > netElementL.sql psql -U postgres -f netElementL.sql eda2.0 }}} Adding a column identifying basin flowing to the sea to rivers {{{#!sql comment ON TABLE portugal.atagua_bacias_bacias_snirh_pc IS 'River basins flowing to the sea'; }}} Trial to join river layer with basin that flows into the sea (this layer "AtAgua_Bacias_bacias_snirh_PC" has the name of basin: "nome") {{{#!sql SELECT rivers.gid, nome FROM portugal.atagua_bacias_bacias_snirh_pc bs -- bassins sea JOIN portugal.rivers ON st_intersects(bs.geom, rivers.geom); SELECT * FROM portugal.rivers LIMIT 10; -- Just a view of the rivers }}} Adding the basin name in the river table: update the table rivers to fill in the column nome using the subquery {{{#!sql ALTER TABLE portugal.rivers ADD COLUMN nome CHARACTER VARYING(50); UPDATE portugal.rivers SET nome = sub.nome FROM( SELECT rivers.gid, bs.nome FROM portugal.atagua_bacias_bacias_snirh_pc bs -- bassins sea JOIN portugal.rivers ON st_intersects( bs.geom, rivers.geom ) ) sub WHERE rivers.gid = sub.gid;-- 75386 SELECT * FROM portugal.rivers WHERE nome IS NOT NULL; }}} Using pgrouting: createTopology also creates vertices_pgr (target is the nearest point to the sea) {{{#!sql SET search_path TO portugal, public; ALTER TABLE portugal.rivers DROP COLUMN IF EXISTS SOURCE; ALTER TABLE portugal.rivers DROP COLUMN IF EXISTS target; ALTER TABLE portugal.rivers ADD COLUMN SOURCE INTEGER; ALTER TABLE portugal.rivers ADD COLUMN target INTEGER; CREATE INDEX rivers_source_idx ON portugal.rivers(SOURCE); CREATE INDEX rivers_target_idx ON portugal.rivers(target); SELECT pgr_createTopology('rivers', 0.0001, 'geom', 'gid'); -- notes: fails with id_localid if character varying so must be gid }}} Create q view of coastal waters: first it is projected to the portugese projection 3763, and second we just need one big water mass for st_distance to the sea hence the UNION {{{#!sql DROP VIEW european_wise2008.coastal_waters_3763; CREATE OR replace VIEW european_wise2008.coastal_waters_3763 AS SELECT st_union( st_transform( the_geom, 3763 ) ) AS geom FROM european_wise2008.costal_waters WHERE cty_id = 'PT'; }}} Creating a table of downstream points {{{#!sql DROP TABLE IF EXISTS portugal.downstream_points; CREATE TABLE portugal.downstream_points AS (SELECT target, nome, TRUE AS at_sea, v.the_geom FROM(SELECT * FROM rivers) r JOIN portugal.rivers_vertices_pgr v ON v.id = r.target WHERE nextdownid = 0);--466 ALTER TABLE portugal.downstream_points ADD CONSTRAINT c_pk_target PRIMARY KEY(target); ALTER TABLE portugal.downstream_points ADD COLUMN distance NUMERIC; -- ALTER TABLE portugal.downstream_points ADD INDEX -- DROP INDEX portugal.downsteam_points_thegeom_idx; CREATE INDEX downstream_points_thegeom_idx ON portugal.downstream_points USING gist(the_geom); }}} Some of the basins are endoreic (red points), not flowing to the sea [[Image(source:eda/data/Docs/trac/sudoang/downstream_points.png, 200px)]][[BR]] So, we calculate distance to the sea and only choose distance <500m as AT_SEA = TRUE {{{#!sql DROP TABLE IF EXISTS portugal.coastalwater_union; CREATE TABLE portugal.coastalwater_union AS SELECT st_union(geom) FROM european_wise2008.coastal_waters_3763; UPDATE portugal.downstream_points SET distance = st_distance(cw.geom, the_geom) FROM european_wise2008.coastal_waters_3763 cw; UPDATE portugal.downstream_points SET at_sea = FALSE WHERE distance > 500; }}} We start with "Lis e Ribeiras Costeiras" as an example, where there are 26 downstream segments {{{#!sql SELECT * FROM downstream_points WHERE nome = 'Lis e Ribeiras Costeiras' --26 }}} ||target||nome||at_sea||the_geom||distance|| ||42602||Lis e Ribeiras Costeiras||f||0101000020B30E00008007CE19EFD0E8C000A8A44E10CDBAC0||29907.6770989189|| ||39176||Lis e Ribeiras Costeiras||t||0101000020B30E0000C0DCB5845C29F0C000470378AB0CE340||0|| ||39812||Lis e Ribeiras Costeiras||t||0101000020B30E000000FBCBEE7D90F0C0005A643BFFC8E040||15.2036257782544|| ||39365||Lis e Ribeiras Costeiras||t||0101000020B30E0000C0627FD93C3FF0C0006ADE71AC8FE240||2.44673969528912|| ||39433||Lis e Ribeiras Costeiras||t||0101000020B30E00004050FC18DD4EF0C000E25817ED31E240||0|| ||39649||Lis e Ribeiras Costeiras||t||0101000020B30E0000802BF6975D74F0C00063EE5AEE69E140||0|| ||39815||Lis e Ribeiras Costeiras||t||0101000020B30E000040A835CD8D85F0C000D3BCE30E09E140||0|| ||39601||Lis e Ribeiras Costeiras||t||0101000020B30E0000C08BDB68BD64F0C000F697DD0DB8E140||0|| ||40528||Lis e Ribeiras Costeiras||t||0101000020B30E000040431CEB6EFFF0C0009ECDAA63EBDC40||9.27728846536787|| ||40005||Lis e Ribeiras Costeiras||t||0101000020B30E0000803255309EACF0C00022FDF6CF3AE040||0|| ||40154||Lis e Ribeiras Costeiras||t||0101000020B30E0000C0F5285CAEBAF0C00010C7BA60C0DF40||10.7529801188864|| ||40328||Lis e Ribeiras Costeiras||t||0101000020B30E00008009F9A0AED3F0C000941804E2A0DE40||0|| ||40371||Lis e Ribeiras Costeiras||t||0101000020B30E00000068226CEEC0F0C0004CC807A17BDF40||2.92664241883266|| ||40616||Lis e Ribeiras Costeiras||t||0101000020B30E000040D8F0F41E04F1C000D2DEE043B6DC40||8.5398098285145|| ||73504||Lis e Ribeiras Costeiras||t||0101000020B30E000000BF0E9CF396EEC0001F85EB7D36E840||0.715212896922981|| ||73603||Lis e Ribeiras Costeiras||t||0101000020B30E00000073D71234B6EEC000A3923A6FC4E740||16.3837915282055|| ||73690||Lis e Ribeiras Costeiras||t||0101000020B30E000080C2F52874BCEEC0006C787A8FAEE740||19.2806257837849|| ||73814||Lis e Ribeiras Costeiras||t||0101000020B30E0000809DEFA7D4DEEEC000EFC9C3F03AE740||66.9101252528966|| ||73894||Lis e Ribeiras Costeiras||t||0101000020B30E000000DFE00B15FEEEC0005F07CE31DDE640||25.1780872038218|| ||74140||Lis e Ribeiras Costeiras||t||0101000020B30E000080499D807684EFC000151DC9056EE540||0|| ||74397||Lis e Ribeiras Costeiras||t||0101000020B30E000000508D979687EFC0006DE7FB455BE540||16.3373603370103|| ||74423||Lis e Ribeiras Costeiras||t||0101000020B30E00008087855AD7BFEFC000CCEEC927A9E440||0|| ||74703||Lis e Ribeiras Costeiras||t||0101000020B30E0000C0EEC9C3FEE7F0C0003AB4C8E2D8DD40||4.00492099466328|| ||74866||Lis e Ribeiras Costeiras||t||0101000020B30E0000C0D32B658F34F1C0002C6519E66DDA40||4.5315444659325|| ||75021||Lis e Ribeiras Costeiras||t||0101000020B30E0000C04B3789BF45F1C00016FBCBA6B8D940||8.50403554270985|| ||75270||Lis e Ribeiras Costeiras||t||0101000020B30E000000E78C289886F1C00070F085892BD740||8.23462515206324|| The target 75270 is the largest flow of this basin Be careful! The target and source attributes are calculating randomly (so this example can be a different number for you) [[Image(source:eda/data/Docs/trac/sudoang/lis_ribeiras_example.png, 300px)]][[BR]] We select everything in basin 75270 The following function (recusrsive) selects everyting in the basin by a recursive query, using source and target created by the pg_routing algorythm (it is much quicker than a spatial recursive) {{{#!sql DROP function if exists portugal.upstream_segments(INTEGER); CREATE FUNCTION portugal.upstream_segments(INTEGER) RETURNS TABLE(source integer, target integer) AS $$ BEGIN RETURN QUERY WITH RECURSIVE parcours AS( SELECT r0.source, r0.target FROM portugal.rivers r0 WHERE r0.target = $1 UNION SELECT r1.source, r1.target FROM parcours JOIN portugal.rivers r1 ON parcours.source = r1.target) SELECT * FROM parcours; END $$ LANGUAGE plpgsql VOLATILE; -- select * from portugal.upstream_segments(75270); --1s and 733rows of source and target }}} Creating indexes to do the routing function and join the result with the river table to get the id of the river segment {{{#!sql ALTER TABLE portugal.rivers ADD COLUMN chemin ltree; CREATE INDEX chemin_gist_rivers_idx ON portugal.rivers USING GIST(chemin); CREATE INDEX chemin_rivers_idx ON portugal.rivers USING btree(chemin); DROP function if exists portugal.get_path(integer,integer); CREATE OR REPLACE FUNCTION portugal.get_path(_from integer, _to integer) RETURNS SETOF ltree AS $$ BEGIN RETURN QUERY SELECT text2ltree(string_agg(gid::text, '.')) AS gid_list FROM (SELECT gid FROM pgr_dijkstra('SELECT gid as id, source, target, st_length(geom) as cost FROM portugal.rivers', _from, _to, FALSE) pt JOIN portugal.rivers r ON pt.edge = r.gid order by seq) sub; END $$ LANGUAGE plpgsql VOLATILE; }}} Different options to calculate topology {{{#!sql select portugal.get_path(75270, source) from portugal.upstream_segments(75270); --5:09 select portugal.get_path(57647, source) from portugal.rivers where nome = 'Lis e Ribeiras Costeiras'; --5:59 UPDATE portugal.rivers set chemin=get_path(75270,u.source) from portugal.upstream_segments(75270) u where rivers.source = u.source; -- 5:13 }}} Pass a vector and insert the values: -- this function writes the path for every basin in an aera (set by the name of the hydrographic region, nome) [[BR]] -- for this, we need to pass the function get_path(target, source) where the target is the node at the estuary [[BR]] -- and the source is all the nodes from the basin connected to that basin [[BR]] -- To achieve the connection to all elements in the basin we used the recursive function upstream_segments [[BR]] Details : -- this function uses a cursor, we could have used a for loop for the cursor [[BR]] -- here following an example http://www.postgresqltutorial.com/plpgsql-cursor/ [[BR]] -- we open the cursor, and loop throught it, and exist at the end with EXIT WHEN NOT FOUND [[BR]] -- We have struggled because we forget to retreive the columns using the * so no column was returned [[BR]] {{{#!sql DROP function if exists portugal.write_chemin(TEXT); CREATE OR REPLACE FUNCTION portugal.write_chemin(_nome text) RETURNS integer AS $$ DECLARE current_count integer default 0; the_downstream_point RECORD; cur_target integer; cur_down CURSOR(_text text) FOR SELECT * FROM portugal.downstream_points WHERE nome = _nome; BEGIN -- Open the cursor OPEN cur_down(_nome); LOOP -- fetch row one by one into the_downstream_point FETCH cur_down INTO the_downstream_point; -- exit when no more row to fetch EXIT WHEN NOT FOUND; current_count := current_count+1; cur_target := the_downstream_point.target; -- raise notice for now RAISE NOTICE 'target :%',cur_target; -- create the chemin for this target and all upstream segments UPDATE portugal.rivers set chemin=get_path(cur_target,u.source) from portugal.upstream_segments(cur_target) u where rivers.source = u.source; END LOOP; -- Close the cursor CLOSE cur_down; RETURN current_count; END; $$ LANGUAGE plpgsql; select portugal.write_chemin('Lis e Ribeiras Costeiras'); select * from portugal.downstream_points where nome = 'Lis e Ribeiras Costeiras'; }}}