| Version 92 (modified by cedric, 15 years ago) (diff) |
|---|
back to first page ..
back to CookBook Eda
back to Obstacle pressure
Join ROE - CCM ¶
Putting geobs2010 into the CCM ¶
Initially, the database was located in another table, below we just copy the table shema roe into the eda2.0 database
I have changed the name geobs into geobs2010 with the new version of the ROE.
REM pour Cédric CD C:\Documents and Settings\cedric\Mes documents\Migrateur\eda\Barrages C:\"Program Files"\PostgreSQL\8.4\bin\pg_dump -U postgres -p 5433 -n geobs2010 ROE> roe_schema.sql C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -h localhost -p 5433 -U postgres -f roe_schema.sql REM==================================== REM Pour Céline CD D:\Celine Jouanin\ROE C:\"Program Files"\PostgreSQL\8.4\bin\pg_dump -U postgres -n geobs2010 ROE> roe_schema.sql C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -h localhost -p 5432 -U postgres -f roe_schema.sql
We have to reproject our data ¶
The constraint CONSTRAINT enforce_srid_ref_position_locale CHECK (srid(ref_position_locale) = 27572) will make it impossible to change the srid from the roe (is that right ?... yes) We will create a new column a new constraint and fill this column with reprojected data
select AddGeometryColumn('geobs2010', 'obstacle_referentiel','ref_position_etrs89' , 3035,'POINT',2);
-- this function also creates constraint « enforce_srid_ref_position_etrs89
select SRID (ref_position_locale) FROM geobs2010.obstacle_referentiel; -- all lines 27572 (Lambert II)
update geobs2010.obstacle_referentiel set ref_position_etrs89 =ST_Transform(ref_position_locale,3035);
Select only the dams within France ¶
Extract the data within the french boundaries ¶
this is ticket #45
With ArcGis ¶
Use the french limits (see French boundaries) in RGF93
In !QuantumGis, connect the ROE database and add the obstacle_referentiel layer
Clic right on the layer "Sauvegarder comme shapefile" with Fichiers de type : Shapefiles (*.shp), Codage=UTF-8 > ok
and choose "Système de coordonnées de référence" = RGF93/Lambert93 (because the french boundary layer is in this system).
In !ArcGis, add the two layers :
- ROE_EDA.shape
- GEOFLADept_FR_Corse_AV_L93\DEPARTEMENT.shp
Use Analysis Tools > Extract > Clip
Input Features : ROE_EDA
Clip Features : DEPARTEMENT
If you want you can choose a "XY Tolerance", because some points in the ROE are in France but they are at the limit.
(maybe 300meters?).
With QGis ¶
first change projection to ERTRS 89
previously we have done this with a request in postgis, this is how you can use gdal
look at Cookbook gdal =>paragraph "changing the projection of any source"
write
ogr2ogr -f formatdst_datasource_name src_datasource_name -a_srs <actual projection> -t_srs <reprojected data>
write ogr2ogr -f formatdst_datasource_name src_datasource_name -a_srs <actual projection> -t_srs <reprojected data>
Launch in OSGeo4W command shell
CD C:\Documents and Settings\cedric\Mes documents\Migrateur\eda\couches_SIG\france C:\OSGeo4W\bin\ogr2ogr DEPARTEMENT_3035.shp DEPARTEMENT.shp -s_srs "EPSG:2154" -t_srs "EPSG:3035"
The shape file DEPARTEMENT_3035.shp is created with right projection
Connect the ROE database and add this newly projected obstacle_referential layer
Process
Add the file GEOFLADept_FR_Corse_AV_L93\DEPARTEMENT.shp (layer>add a vector layer)
Tools > Geoprocessing tools > Clip
But this does not allow to select shapes within a distance, we put this base in postgis, here we create shape along with an index
REM Cédric ======================= CD C:\Documents and Settings\cedric\Mes documents\Migrateur\eda\couches_SIG\france C:\"Program Files"\PostgreSQL\8.4\bin\shp2pgsql -s 3035 -I DEPARTEMENT_3035 france.departement> departement.sql C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -h localhost -U postgres -p 5433 -c "DROP schema if exists france CASCADE" C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -h localhost -U postgres -p 5433 -c "CREATE schema france " C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -h localhost -U postgres -p 5433 -f departement.sql REM Céline ======================= CD D:\Celine Jouanin\france C:\"Program Files"\PostgreSQL\8.4\bin\shp2pgsql -s 3035 -I DEPARTEMENT_3035 france.departement> departement.sql C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -h localhost -U postgres -p 5432 -c "DROP schema if exists france" C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -h localhost -U postgres -p 5432 -c "CREATE schema france" C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -h localhost -U postgres -p 5432 -f departement.sql
insertion dans la table geometry column
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type") SELECT '', 'france', 'departement', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom) FROM france.departement LIMIT 1;
mmh so finally I would like a new column saying whether or not the dams are in France
we have 43937 lines, after joining with France we should have less
SELECT count(ref_id) from geobs.obstacle_referentiel -- 43937SELECT count(ref_id) from geobs2010.obstacle_referentiel -- 47078
The following does not work and it makes sense, we should only have one geom of france
SELECT count(ref_id) from geobs.obstacle_referentiel o join france.departement f ON ST_DWithin(f.the_geom, o.ref_position_etrs89, 2000); --49468
SELECT count(ref_id) from geobs2010.obstacle_referentiel o join france.departement f ON ST_DWithin(f.the_geom, o.ref_position_etrs89, 2000); --52986
I first aggregate all departements and then join this glued layer with the dam layer, the distance is 2000
SELECT count(ref_id) from geobs2010.obstacle_referentiel o join (SELECT ST_Union(f.the_geom) as singlegeom FROM france.departement As f) as sub ON ST_DWithin(sub.singlegeom, o.ref_position_etrs89, 2000); --47049
Build a view allowing to know the score of the different dams ¶
this is ticket #46
Joining CCM and Dams ¶
-- Création d'un index sur obstacle_referentiel
-- Création de la table correspondance (roe_ccm_300), je rajoute un the_geom a la requete pour pouvoir voir la table dans Qgis
La distance de projection entre la CCM et la base geobs est choisie à : 300 m (comme pour l'extraction des données BDMAP sur la CCM)
Attention à l'ordre des colonnes dans la clause group by, si un predicat min est utilisé il faut qu'il apparaisse en deuxième dans la clause
The table should contain
- gid
- dam_height
- dam_score (for the moment defaut is zero)
CREATE INDEX indexgeobs ON geobs2010.obstacle_referentiel USING GIST (ref_position_etrs89 GIST_GEOMETRY_OPS); DROP TABLE IF EXISTS geobs2010.roe_ccm_300; CREATE TABLE geobs2010.roe_ccm_300 as ( SELECT distinct on (ref_id) ref_id, gid, min(distance) as distance, height, score,nbdams,ref_position_etrs89 as the_geom FROM ( SELECT ref_id, gid ,CAST(distance(r.the_geom, b.ref_position_etrs89) as decimal(15,1)) as distance ,b.ref_position_etrs89, CASE WHEN ref_hauteur_chute>0 then ref_hauteur_chute WHEN ref_hauteur_chute=0 then ref_hauteur_chute ELSE ref_hauteur_terrain END AS height, 0 As score, 1 AS nbdams -- pour jointure ultérieure FROM geobs2010.obstacle_referentiel As b INNER JOIN ccm21.riversegments r ON ST_DWithin(r.the_geom, b.ref_position_etrs89,300) -- WHERE b.goodproj IS TRUE -- Attention c'est faux ORDER BY ref_id) AS sub GROUP BY ref_id, distance, gid,height,score,nbdams ,the_geom ); -- mise à jour de la table geometry_columns select Probe_Geometry_Columns(); -- si la commande ci dessus ne marche pas ... INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type") SELECT '', 'geobs2010', 'roe_ccm_300', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom) FROM geobs2010.roe_ccm_300 LIMIT 1; -- creation d'index, clé primaire, et constraintes qui vont bien alter table geobs2010.roe_ccm_300 add column id serial PRIMARY KEY; alter table geobs2010.roe_ccm_300 add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2); alter table geobs2010.roe_ccm_300 add CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'POINT'::text OR the_geom IS NULL); alter table geobs2010.roe_ccm_300 add CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035); CREATE INDEX indexgeobsroe_ccm_300 ON geobs2010.roe_ccm_300 USING GIST ( the_geom GIST_GEOMETRY_OPS );
The final script is below
ALTER TABLE geobs2010.obstacle_referentiel ADD COLUMN goodproj boolean; UPDATE geobs2010.obstacle_referentiel set goodproj=FALSE ; UPDATE geobs2010.obstacle_referentiel set goodproj=TRUE where ref_id IN ( SELECT ref_id from geobs2010.roe_ccm_300);
We now have to use, and can display the dams left out in qgis
select * from geobs2010.obstacle_referentiel where goodproj=TRUE
Some statistics about the dam projection ¶
We create a table of regions
DROP TABLE if exists france.region; CREATE TABLE france.region as( SELECT ST_Union(the_geom) as the_geom,code_reg FROM france.departement GROUP BY code_reg); SELECT geometrytype(the_geom) from france.region; Alter table france.region ADD column reg_id serial PRIMARY KEY; alter table france.region add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2); alter table france.region add CONSTRAINT enforce_geomtype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR geometrytype(the_geom) = 'POLYGON'::text OR the_geom IS NULL); alter table france.region add CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
Calcul du pourcentage de barrages projetés par région.
DROP TABLE if exists geobs2010.pourcentageprojette; CREATE TABLE geobs2010.pourcentageprojette AS( SELECT 100*JOINT/(JOINT +DISJOINT) as pourcentage, nom_region,the_geom FROM( SELECT nom_region, sum(cast(cast(goodproj AS integer)as NUMERIC)) as joint, sum(cast(cast(NOT goodproj AS integer) AS NUMERIC)) as disjoint, the_geom from france.region f, geobs2010.obstacle_referentiel g WHERE ST_Intersects(f.the_geom,g.ref_position_etrs89) group by nom_region, the_geom)as sub); SELECT geometrytype(the_geom) FROM geobs2010.pourcentageprojette; ALTER TABLE geobs2010.pourcentageprojette ADD column pourc_id serial PRIMARY KEY; alter table geobs2010.pourcentageprojette add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2); alter table geobs2010.pourcentageprojette add CONSTRAINT enforce_geomtype_the_geom CHECK(geometrytype(the_geom) = 'MULTIPOLYGON'::text OR geometrytype(the_geom) = 'POLYGON'::text OR the_geom IS NULL);; alter table geobs2010.pourcentageprojette add CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
Pourcentage de barrages projettés avec un buffer de 300 m
(Avec un buffer de 1000, 2000m j'ai toujours les mêmes pourcentages pour chaque région, soit c'est normal soit la requête n'a pas fonctionné).
When we compare the geobs2010.roe_ccm_300 and the geobs2010.obstacle_referentiel :
in red the dams in geobs2010.obstacle_referentiel which are not well projected
in violet the dams which are "well" projected and projected in the CCM layer
With geobs2010.roe_ccm_300, there are two dams (in violet in the East) which are projected in the CCM layer but aren't in France :
ref_id =21164 an ref_id=49967
the table is used in the loaddb method of BaseEdaCCMriversegmentsROE
SELECT sum(height) AS c_height,sum(score) AS c_score,sum(nbdams) AS nbdams,r.gid AS gid FROM geobs2010.roe_ccm_300 j RIGHT JOIN ccm21.riversegments r ON r.gid=j.gid where wso_id in (select wso_id from france.wso where area='France') GROUP BY r.gid ORDER BY r.gid ;
Cumulated number of dams
Cumulated height of dams
