back to start page [[..]]

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

{{{#!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 ==

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

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

{{{#!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 ==
{{{#!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 ==

{{{#!sql
-- 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

-- Bacias dos troços de linha de água GeoCodificadas (basinorder corresponds to idlocalid from rede hordrografica)
shp2pgsql -s 3763 -I vw_baccod_25k_ptcont.shp portugal.vw_baccod_25k_ptcont > vw_baccod_25k_ptcont.sql
psql -U postgres -f vw_baccod_25k_ptcont.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


-- Creating a table corresponding to MINHO
set search_path to portugal ;

-- Table of basins
drop table if exists minho_cuenca;
create table minho_cuenca as select * from vw_baccod_25k_ptcont where basinorder like '181%';

CREATE INDEX minho_cuenca_geom_idx
  ON minho_cuenca
  USING gist
  (geom);

-- Table of rivers
drop table if exists minho_rio;
create table minho_rio as select * from netElementL where id_localid like '181%';

CREATE INDEX minho_rio_geom_idx
  ON minho_rio
  USING gist
  (geom);

-- Configure the table of rivers to use pgrouting: creating source and target
ALTER TABLE minho_rio DROP COLUMN source;
ALTER TABLE minho_rio DROP COLUMN target;
ALTER TABLE minho_rio ADD COLUMN source integer;
ALTER TABLE minho_rio ADD COLUMN target integer;
CREATE INDEX minho_rio_source_idx ON minho_rio (source);
CREATE INDEX minho_rio_target_idx ON minho_rio (target);

---------------------------------------------------------------------
-- Using pgrouting: createTopology also creates vertices_pgr
------------------------------------------------------------------
set search_path to portugal, public;
SELECT pgr_createTopology('minho_rio', 0.0001, 'geom', 'gid');

SELECT  * from minho_rio limit 10

}}}