Noeud - parcours RHT – EDA

Context Navigation

Version 33 (modified by celine, 13 years ago) (diff)
--

back to first page..
back to RHT

Nouvelle version d'octobre 2011 du RHT et topologie (fnode,tnode)

select count(*) from rht.rht ---114601  (version1 avec pb de topologie);
select count(*) from rht.rhtvs2  ---114564 (version octobre);

select id_drain from rht.rht where id_drain not in (select id_drain from rht.rhtvs2);  ---37 id_drain supprimés
select id_drain from rht.rhtvs2 where id_drain not in (select id_drain from rht.rht_topology);  --- 0 id_drain crées

Nouvelle version d'octobre 2011 avec suppression des problèmes de topologie rencontrés (Loire, Seine Normandie) et amélioration du réseau.

INTEGRATION DU RHT et sa TOPOLOGIE

Les données ont été préalablement projetées en 3035 sous ArcGis
Intégration de la table avec id_drain, fnode, tnode

D:
cd D:\CelineJouanin\RHT_Estimkart\RHT_October_newversion
C:\"Program Files"\PostgreSQL\8.4\bin\shp2pgsql -s 3035 -c -g the_geom -W LATIN1 -I rhtvs2_3035.shp rhtvs2 > rhtvs2.sql 
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0_RHT -h localhost -U postgres -p 5432 -f rhtvs2.sql 

ALTER TABLE rhtvs2 SET SCHEMA rht;

-- Rename column
ALTER TABLE rht.rhtvs2 RENAME COLUMN fnode_ TO fnode;
ALTER TABLE rht.rhtvs2 RENAME COLUMN tnode_ TO tnode;
ALTER TABLE rht.rhtvs2 RENAME COLUMN lg TO length;

-- Change type
alter table rht.rhtvs2 alter column id_drain type integer;
alter table rht.rhtvs2 alter column fnode type integer;
alter table rht.rhtvs2 alter column tnode type integer;

-- Constraints
ALTER TABLE rht.rhtvs2 DROP CONSTRAINT rhtvs2_pkey;
ALTER TABLE rht.rhtvs2 ADD CONSTRAINT  pkey_id_drain PRIMARY KEY (id_drain);

-- Create an index on rht
CREATE INDEX indexrhtvs2 ON rht.rhtvs2
  USING GIST ( the_geom GIST_GEOMETRY_OPS );

DROP INDEX IF EXISTS rht.indexrhtvs2;
CREATE INDEX indexrhtv2
ON rht.rhtvs2
USING btree (id_drain);

Noeuds mer RHT

drop table if exists rht.noeudmervs2;
create table rht.noeudmervs2 as select * from rht.rhtvs2 where tnode not in (select fnode from rht.rhtvs2);
CREATE INDEX indexmer ON rht.noeudmervs2 USING btree (id_drain);
CREATE INDEX indexmergeom ON rht.noeudmervs2 USING GIST ( the_geom GIST_GEOMETRY_OPS );
COMMENT ON table rht.noeudmervs2 IS 'Table noeudmervs2 to know the downstream node from the sea';

select count(*) from rht.noeudmervs2; ---1082 noeudmer  (version 1: 1114)

Intégration des noeuds mer dans la table rhtvs2.

alter table rht.rhtvs2 add column noeudmer boolean default FALSE;
update rht.rhtvs2 set noeudmer=TRUE where id_drain in (select id_drain from rht.noeudmervs2); --1082 noeudmer

Nextdownid

create table rht.fromnodetonodevs2 as (select taval.id_drain as id_drainaval, tamont.id_drain as id_drainamont from rht.rhtvs2 tamont join rht.rhtvs2 taval on tamont.tnode=taval.fnode);
alter table rht.rhtvs2 add column nextdownid integer;
update rht.rhtvs2 SET nextdownid=id_drainaval from rht.fromnodetonodevs2 ft where id_drain=ft.id_drainamont;

select  * from rht.rhtvs2 where id_drain=nextdownid --0lignes
update rht.rhtvs2 set nextdownid=-9999 where id_drain=nextdownid --0 lignes

select * from rht.rhtvs2 where noeudmer='t' ---nextdownid null
update rht.rhtvs2 set nextdownid=-9999 where noeudmer='t' --1082 lignes

Noeud source

alter table rht.rhtvs2 add column noeudsource boolean default FALSE;
-- identification des noeuds source
update rht.rhtvs2 set noeudsource=TRUE where id_drain in (
select id_drain from rht.rhtvs2 except (
select rt2.id_drain from rht.rhtvs2 rt1 join rht.rhtvs2 rt2 on rt1.nextdownid=rt2.id_drain)
);-- 57493 lignes

Parcours Distance_sea

Le parcours est fait dans R par la méthode distance_sea.R voir main_RHT_France.R

--pour des raisons d'efficacité on sauve un texte depuis R et on recharge en sql
--Sauvegarde sous R
--rht<-distance_sea(rht)
--write.table(rht@data[,c("id_drain","ltree_from_sea","dmer")],file=str_c(datawd,"/dataEDArht/rht_ltree_distancevs.csv"),sep=";",row.names=FALSE)
--Restauration en sql
create table rht.temp_ltree_distance(
        id_drain integer primary key,
        ltreetext text,
        distance bigint);
alter table rht.temp_ltree_distance alter column distance type double precision;
copy rht.temp_ltree_distance from 'E:/workspace/EDAdata/dataEDArht/rht_ltree_distancevs2.csv' with csv header delimiter as ';' NULL AS 'NA'; --114564 lines

-- lancer le script ltree.sql qui se trouve dans postgres/share
alter table rht.rhtvs2 add column chemin ltree;
update rht.rhtvs2 t set chemin=sub.chemin from (
select text2ltree(ltreetext) as chemin,id_drain from rht.temp_ltree_distance) as sub
where t.id_drain=sub.id_drain; --114564 lines

-- mise à jour des distances mer
alter table rht.rhtvs2 add column dmer numeric;
update rht.rhtvs2 t set dmer=sub.distance from (
select distance,id_drain from rht.temp_ltree_distance) as sub
where t.id_drain=sub.id_drain
and dmer is null;  --114564 lignes

Identification des différentes zones, manche, Bretagne, Méditerranée, Golfe de Gascogne, Est France, Pyrénnées

alter table rht.noeudmervs2 add column zonegeo character varying(30);
update rht.noeudmervs2 n set zonegeo=location from(
select location,id_drain from rht.noeudmer_polygon p
join  rht.noeudmervs2 n on st_intersects(p.the_geom,n.the_geom)) s
where s.id_drain=n.id_drain; -- 1082
update rht.noeudmervs2 set zonegeo='Golfe de Gascogne' where zonegeo='Biscaye';  --188 lignes
select * from rht.noeudmervs2 where zonegeo is null; -- zero lignes

Joining ROE with RHT

Projection spatiale des obstacles (points) sur le RHT (lignes) avec un buffer de 300 m

DROP TABLE IF EXISTS rht.rhtvs2_roev2;
CREATE TABLE rht.rhtvs2_roev2 as (
        SELECT distinct on (id_roe) id_roe, id_drain, min(distance) as distance, the_geom FROM (
               SELECT id_roe, id_drain ,CAST(distance(r.the_geom, s.the_geom) as  decimal(15,1)) as distance, s.the_geom 
               FROM roe_v2.roe_v2 As s
               INNER JOIN  rht.rhtvs2 r ON ST_DWithin(r.the_geom, s.the_geom,300)
               WHERE s.the_geom IS NOT NULL
               ORDER BY id_roe) AS sub 
        GROUP BY id_roe, distance,id_drain, the_geom  ---49269 lines
); 
alter table rht.rhtvs2_roev2 add column id serial;
-- mise à jour de la table geometry_columns
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'rht', 'rhtvs2_roev2', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM rht.rhtvs2_roev2 LIMIT 1;

-- creation d'index, clé primaire, et constraintes qui vont bien
alter table rht.rhtvs2_roev2 add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table rht.rhtvs2_roev2 add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'POINT'::text OR the_geom IS NULL);
alter table rht.rhtvs2_roev2 add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
alter table rht.rhtvs2_roev2 ADD CONSTRAINT pk_idroe PRIMARY KEY(id);
alter table rht.rhtvs2_roev2 drop CONSTRAINT pk_idroe;
CREATE INDEX indexroev2_rhtvs2 ON rht.rhtvs2_roev2
  USING GIST ( the_geom GIST_GEOMETRY_OPS );

Joining BDMAP with RHT

Projection spatiale des points BDMAP sur la couche station_geography (bdmap2009).

-- creation de la table bdmap2009.bdmap_rhtvs2, je rajoute un the_geom a la requete pour pouvoir voir la table dans Qgis  
DROP TABLE IF EXISTS bdmap2009.bdmap_rhtvs2;
CREATE TABLE bdmap2009.bdmap_rhtvs2 as (
        SELECT distinct on (st_codecsp) st_codecsp, id_drain, min(distance) as distance, the_geom FROM (
               SELECT st_codecsp, id_drain ,CAST(distance(r.the_geom, s.the_geom) as  decimal(15,1)) as distance,s.the_geom 
               FROM bdmap2009.stationsp2 As s
               INNER JOIN  rht.rhtvs2 r ON ST_DWithin(r.the_geom, s.the_geom,300)
               WHERE s.the_geom IS NOT NULL
               ORDER BY st_codecsp) AS sub 
        GROUP BY st_codecsp, distance,id_drain, the_geom
);
alter table bdmap2009.bdmap_rhtvs2 add column id serial;
-- mise à jour de la table geometry_columns
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'bdmap2009', 'bdmap_rhtvs2', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM bdmap2009.bdmap_rhtvs2 LIMIT 1;

-- creation d'index, clé primaire, et constraintes qui vont bien
alter table bdmap2009.bdmap_rhtvs2 add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table bdmap2009.bdmap_rhtvs2 add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'POINT'::text OR the_geom IS NULL);
alter table bdmap2009.bdmap_rhtvs2 add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
alter table bdmap2009.bdmap_rhtvs2 ADD CONSTRAINT pk_id PRIMARY KEY(id);
CREATE INDEX indexbdmap_rht ON bdmap2009.bdmap_rhtvs2
  USING GIST ( the_geom GIST_GEOMETRY_OPS );

#A FINIR ci-dessous !!

Comment aller chercher les infos sur le nombre de barrages et la distance mer, les données sont dans la ccm ?

installation de dblink (je suis aller chercher le code sql dans contrib), dblink s'installe dans une base de données PostgreSQL: pgsql/contrib/dblink/dblink.sql

/* cross database query*/
-- je crée une table contenant la ccm.riversegments pour la France avec les colonnes qui m'intéressent.

select dblink_connect('connectionccm','hostaddr=93.20.247.238 port=5432 dbname=eda2 user=postgres password=petromyzon***');
create schema ccm;
drop table if exists ccm.riversegments_france;
create table ccm.riversegments_france as (
SELECT * FROM dblink('connectionccm','SELECT wso1_id,cs_nbdams,cum_len_sea,the_geom  FROM ccm21.riversegments r join europe.wso w on r.wso_id=w.wso_id where area=''France''')
 AS  t(wso1_id int, cs_nbdams integer, cum_len_sea numeric,the_geom geometry)
);

-- je vais chercher dans la ccm les informations sur les distances mer et les cumuls de barrages

select distinct on (zonegeo) zonegeo from rht.noeudmervs2;

alter table rht.noeudmervs2 add column cumnbbar integer;
alter table rht.noeudmervs2 add column dmer numeric;
alter table rht.noeudmervs2 add column wso1_id integer;


update rht.noeudmervs2 set (cumnbbar,dmer,wso1_id)= (jointure.cs_nbdams,jointure.cum_len_sea,jointure.wso1_id) from
        -- resultat final de la jointure spatiale entre ccm et noeudfrontiere
        (select distinct on (id_drain) id_drain, min(distance) as distance, wso1_id,cs_nbdams, cum_len_sea from (
                -- requete intermédiaire dont on prend le min (sub)
                select id_drain, ST_distance(noeudfrontiere.the_geom,ccmrs.the_geom) as distance,ccmrs.*
                FROM ccm.riversegments_france ccmrs join
                (select id_drain, the_geom from rht.noeudmervs2 where noeudmer and zonegeo in ('Est','pyrénées')) as noeudfrontiere
                on ST_DWithin(noeudfrontiere.the_geom,ccmrs.the_geom, 500)) as sub
        group by id_drain, distance, wso1_id,cs_nbdams, cum_len_sea) as jointure
where noeudmer.id_drain=jointure.id_drain; --100 lignes modifiées

=> sauvegarde pour Céline

e: cd E:\IAV\eda\rht pg_dump -U postgres -p 5432 -t rht.rht_topology -t rht.noeudmervs2 -t rht.noeudmervs2_polygon eda2.0_RHT> sauve3tablesnoeudsmer.sql

---Celine d: cd D:\CelineJouanin?\EDA20RHT C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0_RHT -h localhost -U postgres -p 5432 -f sauve3tablesnoeudsmer.sql

Mise à jour de la table rht_topology

alter table rht.rht_topology add column cumnbbar integer;
alter table rht.rht_topology add column dmer numeric;
alter table rht.rht_topology drop column wso1_id ;

update  rht.noeudmervs2 set dmer=0 where dmer is null and noeudmer; -- 983 -- a modifier pour prendre peut être la moitié de la distance du tronçon... a voir avec Céline
select * from rht.noeudmervs2 where dmer>0; -- il n'a pas trouvé de bassin dans les pyrénées... bizarre.... Mais on vivra bien avec !!!
-- mise à jour des données dans rht_topology
update rht.rht_topology set (dmer, cumnbbar)=(sub.dmer, sub.cumnbbar)
from (select id_drain,dmer,cumnbbar from rht.noeudmervs2 where dmer is not null and noeudmer) sub
where sub.id_drain=rht_topology.id_drain; --1083

Distance source (Parcours source)

-- script de calcul des distances sources
set search_path to rht,public;
select cast(ltree2text(subltree(chemin,0,1)) as integer) from rht_topology limit 10
/*
AJOUT D'UN INDEX
ltree supports several types of indexes that can speed up the indicated operators:
    B-tree index over ltree: <, <=, =, >=, >
    GiST index over ltree: <, <=, =, >=, >, @>, <@, @, ~, ? 
*/
CREATE INDEX chemin_btee_index ON rht.rht_topology USING btree(chemin); --27688 ms
VACUUM ANALYSE rht.rht_topology;
CREATE INDEX chemin_gist_index ON rht.rht_topology USING GIST (chemin); -- fait planter mon serveur....
SELECT * from rht.rht_topology where chemin <@ '212340'; -- tronçon source de la Vilaine -- 2360
SELECT chemin from rht.rht_topology where chemin <@ '212340' limit 1
-- La distance source est la distance mer maximale de tous les chemins possibles en ne prenant que le dernier tronçon du chemin.
Select max(dmer) from rht.rht_topology where id_drain in (
Select cast(ltree2text(subpath(chemin, -1, 1))as integer) as chenin_id_drain from rht.rht_topology where chemin ~ '*.212340.*');
--  travail sur les noeuds source et le recalcul des distances mer et source (voir fichier doc joint)
-- en pratique les distances sont calculées a partir du milieu de chaque tronçon, elles prennent par contre toute la distance des tronçons autre que le tronçon considéré.

alter table rht.rht_topology add column noeudsource boolean default FALSE;
-- identification des noeuds source
update rht.rht_topology set noeudsource=TRUE where id_drain in (
select id_drain from rht.rht_topology except (
select rt2.id_drain from rht.rht_topology rt1 join rht.rht_topology rt2 on rt1.nextdownid=rt2.id_drain)
);-- 57533
update rht.rht_topology set dmer=dmer+length/2;-- 114600
alter table rht.rht_topology add column dsource numeric;
update rht.rht_topology set dsource=length/2 where noeudsource;-- 57533
-- voir schéma pour explications


DROP FUNCTION IF EXISTS rht.distance_source();

CREATE OR REPLACE FUNCTION rht.distance_source(min integer, max integer) RETURNS integer AS $$
DECLARE
    ligne integer;
    nbligne integer; 
    id_drains record;
BEGIN
        ligne:=1;
        DROP TABLE IF EXISTS sousselection;
        DROP sequence IF EXISTS  seq;
        CREATE sequence seq;
        EXECUTE 'create temporary table sousselection as SELECT id_drain, nextval(''seq'') as num FROM (SELECT * FROM rht.rht_topology WHERE noeudsource=FALSE ORDER BY id_drain) as sub WHERE id_drain >='|| min || ' AND id_drain<='||max ;
        SELECT COUNT(*) FROM sousselection INTO nbligne;
        WHILE ligne <=nbligne LOOP
                  -- pour vérification décommenter ensuite
                  SELECT * from sousselection where num=ligne INTO id_drains;
                  RAISE NOTICE 'distance_source for %', quote_ident(CAST(id_drains.id_drain as text));
                  EXECUTE 'INSERT INTO rht.distance_source SELECT '||id_drains.id_drain||', MAX(dmer+length/2) as dmer FROM rht.rht_topology WHERE id_drain IN (
                        SELECT CAST(ltree2text(subpath(chemin, -1, 1))as integer) FROM rht.rht_topology where chemin ~  ''*.'||
                        id_drains.id_drain
                        || '.*'' )';  
                ligne=ligne+1;      
        END LOOP;
    RETURN ligne;
END;
$$ LANGUAGE plpgsql;
-- lancement de la requête
DROP TABLE IF EXISTS rht.distance_source;
CREATE TABLE rht.distance_source (id_drain integer,distance  numeric);
select rht.distance_source(1,50000);

-- pour essai pour voir ce qui se passe dans la requête

SELECT 212340,  MAX(dmer+length/2) as dmer FROM rht.rht_topology WHERE id_drain IN (
                        SELECT CAST(ltree2text(subpath(chemin, -1, 1))as integer) FROM rht.rht_topology where chemin ~  '*.212340.*' )
-- il faudra enlever la dmer de 212340

-- POUR LES LANCEMENTS FINAUX IL A FALLU LANCER UN INDEX ET RAJOUTER UNE CLE PRIMAIRE
Alter table rht.rht_topology add constraint c_pk_id_drain primary key (id_drain);
CREATE INDEX id_drain_index  ON rht.rht_topology
  USING btree
  (id_drain);

select rht.distance_source(100000,200000);
select rht.distance_source(200001,300000);
select rht.distance_source(300001,400000);
select rht.distance_source(400001,600000);

-- pour essai pour voir ce qui se passe dans la requête

SELECT 212340,  MAX(dmer+length/2) as dmer FROM rht.rht_topology WHERE id_drain IN (
                        SELECT CAST(ltree2text(subpath(chemin, -1, 1))as integer) FROM rht.rht_topology where chemin ~  '*.212340.*' )


select * from rht.rht_topology t join rht.distance_source s on s.id_drain=t.id_drain where distance is null; --323

-- il faudra enlever la dmer de 212340

-- requête finale pour calculer les distances source (à lancer)
alter table rht.distance_source  add column distance_source numeric;
select count(*) from rht.distance_source; --57067 que les noeuds non source

-- mise à jour des distances source dans la table distance source (il faut encore enlever dmer)
update rht.distance_source d set distance_source = sub.distance_source from (
select t.id_drain, distance-dmer as distance_source from rht.distance_source s join rht.rht_topology t on s.id_drain=t.id_drain) as sub
where sub.id_drain=d.id_drain; --57067

-- copie des valeurs de distance_source dans dsource de la table rht.rht_topology
update rht.rht_topology t set dsource=distance_source
from rht.distance_source d where d.id_drain=t.id_drain; --57067

-- mise à jour des distances_source (pour l'instant valeurs = null) pour les noeuds source (pour les derniers c'est la moitié de la longueur du segment)
update rht.rht_topology t set dsource=length/2 where noeudsource;
-- nombre de 
select count(*)as count, noeudsource from rht.rht_topology where dsource is null group by noeudsource;

Utilisation des donnees de parcours

En pratique nous n'aurons probablement besoin que du noeud le plus amont des parcours comme programmé ci dessus. Mais les scripts ci-dessous pourraient s'avérer utiles.
Ci dessous la fonction utilise les expressions régulières ltree => vecteur => table

/*
regexp_split_to_array(string text, pattern text [, flags text ])        regexp_split_to_array('hello world', E'\\s+')   {hello,world}
SELECT regexp_split_to_array('hello world', E'\\s+');
SELECT regexp_split_to_array('hello.world', E'\\.+');
*/
SELECT regexp_split_to_array(ltree2text(chemin),E'\\.+') from rht.rht_topology where chemin <@ '212340';

/*
http://www.pgsql.cz/index.php/PostgreSQL_SQL_Tricks#Array_to_table
SELECT (ARRAY[1,3,4,6,7])[idx.n] FROM generate_series(1,7) idx(n);
CREATE OR REPLACE FUNCTION unpack(anyarray)
RETURNS SETOF anyelement AS $$ 
SELECT $1[i] 
   FROM generate_series(array_lower($1,1), 
                        array_upper($1,1)) g(i);
$$ LANGUAGE sql STRICT IMMUTABLE;
*/
CREATE OR REPLACE FUNCTION unpack(anyarray)
RETURNS SETOF anyelement AS $$ 
SELECT $1[i] 
   FROM generate_series(array_lower($1,1), 
                        array_upper($1,1)) g(i);
$$ LANGUAGE sql STRICT IMMUTABLE;
-- pour convertir un vecteur en table
select unpack(vecteurchemin) from 
        (SELECT regexp_split_to_array(ltree2text(chemin),E'\\.+') as vecteurchemin from rht.rht_topology where chemin <@ '212340' limit 1) as sub
-- j'avais pas vu il y a déjà une fonction magique unest...
select cast(unnest(vecteurchemin) as integer) as chemin_id_drain from
        (SELECT regexp_split_to_array(ltree2text(chemin),E'\\.+') as vecteurchemin from rht.rht_topology where chemin <@ '212340' limit 1) as sub
-- Fonction qui renvoit les identifiants uniques des tronçons en amont d'un ouvrage


CREATE OR REPLACE FUNCTION troncons(integer)
RETURNS SETOF integer AS $$
select cast(unnest(vecteurchemin) as integer) as chemin_id_drain from
        (SELECT regexp_split_to_array(ltree2text(chemin),E'\\.+') as vecteurchemin from rht.rht_topology where chemin <@ text2ltree(cast ($1 as text))  limit 1) as sub;
$$ LANGUAGE sql STRICT IMMUTABLE;       

SELECT troncons_amont(212340); -- OK

Download in other formats:

Plain Text