wiki:import layers postgis

Version 47 (modified by maria, 7 years ago) (diff)

--

back to start page ..

      1. Tickets
  2. Import layer postgis
    1. first fixing a problem with the database
    2. First trial with two layers
    3. Trials for the pfafsteter system
    4. pg_routing
    5. ltree
    6. Select some bassins to calculate topology
    7. Trial with Portugal layers
    8. Combining postgis and R to calculate topology

Tickets

#90 #97

Import layer postgis

first fixing a problem with the database

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 

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 

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

Our objective is to test the best way to route the river network 

-- 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 

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 

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"""27072707"A"175359"10037294"...

10037 corresponds to the Biscay coast of Spain, 1003714 correspond à Oria 

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 

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 

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);

source:eda/data/Docs/trac/sudoang/oria.png

Creating a table corresponding to UROLA

Table of basins 

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 

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);

source:eda/data/Docs/trac/sudoang/urola.png

Trials for the pfafsteter system

To be able to explore the topology we could use pfafstetter 

/*-----------------------------------------
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'

pfafcuenlngtramo_m
10037146451747
10037146453304
10037146455437

Problem from that step we need to go to another number

10037146451 => 1003714643 

-- 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
where 1037146(47) is immediately upstream from 1037146(457)
source:eda/data/Docs/trac/sudoang/network_topology.PNG

pg_routing

​https://docs.pgrouting.org/2.5/en/pgRouting-concepts.html#description-of-the-edges-sql-query-for-dijkstra-like-functions 

-------------------------------------
-- 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"
157740774"1003714113"
228783340.239109235141783"1003714115"
333916434.71909815275916"1003714117"
4619201381.42299559645920"100371413"
5629561931.41349601216956"10037141511"
6549922751.07786782504992"10037141513"
75310233006.697040024421023"10037141515"
81610283923.516044522871028"10037141517"
9159624550.98432547035962"10037141531"
10478795241.61915518433879"10037141533"
11398805777.24646955956880"10037141535"

ltree

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

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
source:eda/data/Docs/trac/sudoang/test_biscay_coast.png

pgRouting, configuration on table and ltree 

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

In Command prompt(cmd) we need to create a new schema and import the shapefiles to Postgres 

-- 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 

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") 

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 

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) 

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 a 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 

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 

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

source:eda/data/Docs/trac/sudoang/downstream_points.png

So, we calculate distance to the sea and only choose distance <500m as AT_SEA = TRUE 

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 

SELECT * FROM downstream_points WHERE nome = 'Lis e Ribeiras Costeiras' --26
targetnomeat_seathe_geomdistance
42602Lis e Ribeiras Costeirasf0101000020B30E00008007CE19EFD0E8C000A8A44E10CDBAC029907.6770989189
39176Lis e Ribeiras Costeirast0101000020B30E0000C0DCB5845C29F0C000470378AB0CE3400
39812Lis e Ribeiras Costeirast0101000020B30E000000FBCBEE7D90F0C0005A643BFFC8E04015.2036257782544
39365Lis e Ribeiras Costeirast0101000020B30E0000C0627FD93C3FF0C0006ADE71AC8FE2402.44673969528912
39433Lis e Ribeiras Costeirast0101000020B30E00004050FC18DD4EF0C000E25817ED31E2400
39649Lis e Ribeiras Costeirast0101000020B30E0000802BF6975D74F0C00063EE5AEE69E1400
39815Lis e Ribeiras Costeirast0101000020B30E000040A835CD8D85F0C000D3BCE30E09E1400
39601Lis e Ribeiras Costeirast0101000020B30E0000C08BDB68BD64F0C000F697DD0DB8E1400
40528Lis e Ribeiras Costeirast0101000020B30E000040431CEB6EFFF0C0009ECDAA63EBDC409.27728846536787
40005Lis e Ribeiras Costeirast0101000020B30E0000803255309EACF0C00022FDF6CF3AE0400
40154Lis e Ribeiras Costeirast0101000020B30E0000C0F5285CAEBAF0C00010C7BA60C0DF4010.7529801188864
40328Lis e Ribeiras Costeirast0101000020B30E00008009F9A0AED3F0C000941804E2A0DE400
40371Lis e Ribeiras Costeirast0101000020B30E00000068226CEEC0F0C0004CC807A17BDF402.92664241883266
40616Lis e Ribeiras Costeirast0101000020B30E000040D8F0F41E04F1C000D2DEE043B6DC408.5398098285145
73504Lis e Ribeiras Costeirast0101000020B30E000000BF0E9CF396EEC0001F85EB7D36E8400.715212896922981
73603Lis e Ribeiras Costeirast0101000020B30E00000073D71234B6EEC000A3923A6FC4E74016.3837915282055
73690Lis e Ribeiras Costeirast0101000020B30E000080C2F52874BCEEC0006C787A8FAEE74019.2806257837849
73814Lis e Ribeiras Costeirast0101000020B30E0000809DEFA7D4DEEEC000EFC9C3F03AE74066.9101252528966
73894Lis e Ribeiras Costeirast0101000020B30E000000DFE00B15FEEEC0005F07CE31DDE64025.1780872038218
74140Lis e Ribeiras Costeirast0101000020B30E000080499D807684EFC000151DC9056EE5400
74397Lis e Ribeiras Costeirast0101000020B30E000000508D979687EFC0006DE7FB455BE54016.3373603370103
74423Lis e Ribeiras Costeirast0101000020B30E00008087855AD7BFEFC000CCEEC927A9E4400
74703Lis e Ribeiras Costeirast0101000020B30E0000C0EEC9C3FEE7F0C0003AB4C8E2D8DD404.00492099466328
74866Lis e Ribeiras Costeirast0101000020B30E0000C0D32B658F34F1C0002C6519E66DDA404.5315444659325
75021Lis e Ribeiras Costeirast0101000020B30E0000C04B3789BF45F1C00016FBCBA6B8D9408.50403554270985
75270Lis e Ribeiras Costeirast0101000020B30E000000E78C289886F1C00070F085892BD7408.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)

source:eda/data/Docs/trac/sudoang/lis_ribeiras_example.png

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) 

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 

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:

-- For our example "Lis e Ribeiras Costeiras", in the largest flow "75270" 

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)
-- for this, we need to pass the function get_path(target, source) where the target is the node at the estuary
-- and the source is all the nodes from the basin connected to that basin
-- To achieve the connection to all elements in the basin we used the recursive function upstream_segments

Details:

-- this function uses a cursor, we could have used a for loop for the cursor
-- here following an example ​http://www.postgresqltutorial.com/plpgsql-cursor/
-- we open the cursor, and loop throught it, and exist at the end with EXIT WHEN NOT FOUND
-- We have struggled because we forget to retreive the columns using the * so no column was returned 

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;

We run the next function for each portuguese basin (19 basins) to calculate the path of rivers that we will use later in te model 

select portugal.write_chemin('Lis e Ribeiras Costeiras');

Combining postgis and R to calculate topology

After several issues calculating topology in PostGis? (problem with the MB in max_stack_depth and the high calculation time when basins are big), some methods created in R for the CCM layer are used

At first, some parameters are created in sql:

-- nextdownid: estimate the downstream gid located next to each segment
-- dmer: adding zeros for the downstream points (the method will calculate later the distance to the sea) 

-- Summary of tables to see what there are in rivers and downstream_points
SELECT count(*) FROM portugal.rivers;
SELECT * FROM portugal.rivers limit 10;
SELECT count(*) FROM portugal.downstream_points;
SELECT * FROM portugal.downstream_points limit 10;

-- Change the name of "nextdownid" because they use the "hydroid" code and we want to stay with "target"
ALTER TABLE portugal.rivers RENAME COLUMN nextdownid TO nextdownid_hydroid;

-- So, we create a new nextdownid
--ALTER TABLE portugal.rivers DROP COLUMN IF EXISTS nextdownid;
ALTER TABLE portugal.rivers ADD COLUMN nextdownid integer;

UPDATE portugal.rivers SET nextdownid = sub.nextdownid FROM
(SELECT r2.gid AS nextdownid, r1.gid FROM portugal.rivers r1 JOIN portugal.rivers r2 ON r2.source = r1.target 
) sub
WHERE sub.gid = rivers.gid; -- 100765 rows affected, 6.6 secs

CREATE index idx_rivers_nextdownid ON portugal.rivers(nextdownid);


-- Add dmer column to calculate distance to the sea (I do not know if the length_ column from rivers is the same)
ALTER TABLE portugal.rivers ADD COLUMN dmer numeric;
UPDATE portugal.rivers SET dmer = 0 WHERE target IN (SELECT target FROM portugal.downstream_points); --466 rows, 363 msec

This is the portuguese script, but it is the same for Spain (in Spain "nextdownid" did not exist, so it was also created).

Once these parameters are created, we use some methods in R to estimate the distance to the sea and the pathfinding: source:eda/EDA/EDA_Rios/BaseEdaRiosRiversegments.R