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 ccm21_eda database
{{{
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 ccm21_eda -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 geobs ROE> roe_schema.sql
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d ccm21_eda -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 ?)
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 ccm21_eda -h localhost -U postgres -p 5433 -c "DROP schema if exists france CASCADE"
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d ccm21_eda -h localhost -U postgres -p 5433 -c "CREATE schema france "
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d ccm21_eda -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 ccm21_eda -h localhost -U postgres -p 5432 -c "DROP schema if exists france"
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d ccm21_eda -h localhost -U postgres -p 5432 -c "CREATE schema france"
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d ccm21_eda -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
}}}
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
}}}
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
}}}
The final script is below
{{{
ALTER TABLE geobs.obstacle_referentiel ADD COLUMN goodproj boolean;
UPDATE geobs.obstacle_referentiel set goodproj=TRUE where ref_id IN (
SELECT  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));	
}}}
filling the holes
{{{
UPDATE geobs.obstacle_referentiel set goodproj=FALSE where  goodproj is NULL;
}}}
We now have to use, and can display the dams left out in qgis
{{{
select * from geobs.obstacle_referentiel where goodproj=TRUE
}}}
[[Image(source:data/Docs/trac/barrages_cutfrance.png,800px)]]
== 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  riversegments r ON ST_DWithin(r.the_geom, b.ref_position_etrs89,300)
               WHERE b.goodproj IS TRUE
               ORDER BY ref_id) AS sub 
        GROUP BY ref_id, gid, the_geom
); --53s

-- mise à jour de la table geometry_columns
alter table geobs.correspondance add column id serial PRIMARY KEY;
select Probe_Geometry_Columns(); -- 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 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 );
}}}

== 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_Collect(the_geom) as the_geom,code_reg FROM france.departement GROUP BY code_reg);
Alter table france.region ADD column reg_id serial PRIMARY KEY;