back to first page [..] [[BR]]
back to ["CookBook Eda"] [[BR]]


= Join ROE - CCM =
== Putting geobs 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 change 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 -n geobs 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('geobs', 'obstacle_referentiel','ref_position_etrs89' , 3035,'POINT',2);
-- this function also creates constraint « enforce_srid_ref_position_etrs89 
select SRID (ref_position_locale) FROM geobs.obstacle_referentiel; -- all lines 27572 (Lambert II)
update geobs.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[[BR]]
In __!QuantumGis__, connect the ROE database and add the obstacle_referentiel layer[[BR]]
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).[[BR]]
In __!ArcGis__, add the two layers :[[BR]]
- ROE_EDA.shape[[BR]]
- GEOFLADept_FR_Corse_AV_L93\DEPARTEMENT.shp[[BR]]
Use Analysis Tools > Extract > Clip[[BR]]
Input Features : ROE_EDA[[BR]]
Clip Features : DEPARTEMENT[[BR]]

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?).[[BR]]

==== With QGis ====

first change projection to ERTRS 89[[BR]]
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"[[BR]]
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>[[BR]]
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[[BR]]
Connect the ROE database and add this newly projected obstacle_referential layer [[BR]]
Process 
Add the file GEOFLADept_FR_Corse_AV_L93\DEPARTEMENT.shp (layer>add a vector layer) [[BR]]
Tools > Geoprocessing tools > Clip[[BR]]
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[[BR]]
{{{
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 [[BR]]
 we have 43937 lines, after joining with France we should have less
{{{
SELECT count(ref_id) from geobs.obstacle_referentiel -- 43937
}}}
{{{
SELECT 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 geobs.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) --43909
}}}

[[Image(source:data/Docs/trac/barrages_cutfrance.png,600px)]][[BR]]

== 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, je rajoute un the_geom a la requete pour pouvoir voir la table dans Qgis 
{{{
CREATE INDEX indexgeobs ON geobs.obstacle_referentiel  USING GIST (ref_position_etrs89 GIST_GEOMETRY_OPS);
DROP TABLE IF EXISTS geobs.correspondance;
CREATE TABLE geobs.correspondance as (
        SELECT distinct on (ref_id) ref_id, gid, min(distance) as distance, 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
               FROM geobs.obstacle_referentiel As b
               INNER JOIN  ccm21.riversegments r ON ST_DWithin(r.the_geom, b.ref_position_etrs89,1000)
--             WHERE b.goodproj IS TRUE -- Attention c'est faux
               ORDER BY ref_id) AS sub 
        GROUP BY ref_id, gid, the_geom
); --53s

-- 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 '', 'geobs', 'correspondance', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM geobs.correspondance LIMIT 1;

-- creation d'index, clé primaire, et constraintes qui vont bien
alter table geobs.correspondance add column id serial PRIMARY KEY;
alter table geobs.correspondance add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table geobs.correspondance add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'POINT'::text OR the_geom IS NULL);
alter table geobs.correspondance add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);

CREATE INDEX indexgeobscorrespondance ON geobs.correspondance
  USING GIST ( the_geom GIST_GEOMETRY_OPS );
}}}
The final script is below
{{{
ALTER TABLE geobs.obstacle_referentiel ADD COLUMN goodproj boolean;
UPDATE geobs.obstacle_referentiel set goodproj=FALSE ;
UPDATE geobs.obstacle_referentiel set goodproj=TRUE where ref_id IN (
SELECT  ref_id from geobs.correspondance);
}}}
We now have to use, and can display the dams left out in qgis
{{{
select * from geobs.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 geobs.pourcentageprojette;
CREATE TABLE geobs.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, geobs.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 geobs.pourcentageprojette;
ALTER TABLE geobs.pourcentageprojette ADD column pourc_id serial PRIMARY KEY;
alter table geobs.pourcentageprojette  add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table geobs.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 geobs.pourcentageprojette  add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
}}}
Pourcentage de barrages projettés avec un buffer de 300 m

[[Image(source:data/Docs/trac/pourc_projette.jpg)]]

(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é).[[BR]]


When we compare the geobs.correspondance and the geobs.obstacle_referentiel :[[BR]]
[[Image(source:data/Docs/trac/ROE/Obstacle referentiel et correspondance.png,800px)]]
in red the dams in geobs.obstacle_referentiel which are not well projected
in violet the dams which are "well" projected and projected in the CCM layer
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