back to first page ..
back to CookBook Eda

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 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 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 geobs2010.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 geobs2010.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) --43909

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 geobs2010.obstacle_referentiel  USING GIST (ref_position_etrs89 GIST_GEOMETRY_OPS);
DROP TABLE IF EXISTS geobs2010.correspondance;
CREATE TABLE geobs2010.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 geobs2010.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 '', 'geobs2010', 'correspondance', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM geobs2010.correspondance LIMIT 1;

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

CREATE INDEX indexgeobscorrespondance ON geobs2010.correspondance
  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.correspondance);

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