back to first page ..
back to CookBook Eda

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
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 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 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

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 -- 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

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

Build a view allowing to know the score of the different dams

this is ticket #46

Joining CCM and Dams

this is ticket #47