Changes between Version 38 and Version 39 of import layers postgis

Timestamp:

Aug 21, 2018 12:05:26 PM (7 years ago)

Author:

maria

Comment:

--

import layers postgis

@@ -414 +414 @@
 == Trial with Portugal layers ==
 
-{{{#!sql
--- In Command prompt(cmd) we need to create a new schema and import the shapefiles to Postgres
+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
@@ -424 +425 @@
 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
-
-}}}
+}}}
+
+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';
+}}}