back to first page [..][[BR]]
back to Recipes for EDA ["CookBook Eda"]

= How to join the BDMAP and CCM version 2 =
see Joining BDMAP and the CCM layer ["CookBook join BDMAP_CCM"][[BR]] for the first version. [[BR]]
The trouble is that some dams or electrofishing stations are sometimes badly projected [[BR]]
We should start at step 6 [[BR]]

== 1. putting bd_carthage into postgis ==
below an example on how to load from latin1 database to postgis
only hylcov_arc
{{{
cd C:\eda\couches_SIG\BDCarthage09
C:\"Program Files"\PostgreSQL\8.4\bin\psql -p 5433 -U postgres -c "CREATE DATABASE bd_carthage  WITH OWNER = postgres  template=postgis ENCODING = 'UTF8'  LC_COLLATE = 'French_France.1252'  LC_CTYPE = 'French_France.1252' CONNECTION LIMIT = -1;"
C:\"Program Files"\PostgreSQL\8.4\bin\psql -p 5433 -U postgres -c "COMMENT ON DATABASE bd_carthage IS 'database bd_carthage2009" 
C:\"Program Files"\PostgreSQL\8.4\bin\shp2pgsql -s 3035 -c -g the_geom -W LATIN1 -I hylcov_arc.shp hylcov_arc> hylcov_arc.sql 
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d bd_carthage -h localhost -U postgres -p 5433 -f hylcov_arc.sql 
}}}
~~We also need the nodes from bd_carthage~~ use downstream node instead

{{{
cd C:\eda\couches_SIG\BDCarthage09
C:\"Program Files"\PostgreSQL\8.4\bin\shp2pgsql -s 3035 -c -g the_geom -W LATIN1 -I hylcov_node.shp hylcov_node> hylcov_node.sql 
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d bd_carthage -h localhost -U postgres -p 5433 -f hylcov_node.sql 
}}}

changing the projections (if the layer is not already in 3035 the first four lines can be avoided by putting -s 2154 in the previous box

{{{
#!sql
-- les projections n'étaient pas en 3035, erreur ci dessus que je corrige
alter table hylcov_node drop constraint enforce_srid_the_geom ;
alter table hylcov_arc drop constraint enforce_srid_the_geom ;
update hylcov_node set the_geom=ST_SetSRID(the_geom, 2154);
update hylcov_arc set the_geom=ST_SetSRID(the_geom, 2154);
update hylcov_arc set the_geom=ST_Transform(the_geom,3035);
update hylcov_node set the_geom=ST_Transform(the_geom,3035);
alter table hylcov_node add constraint enforce_srid_the_geom CHECK (st_srid(the_geom) = 3035);
alter table hylcov_arc add constraint enforce_srid_the_geom CHECK (st_srid(the_geom) = 3035);
}}}


== 2. getting distance sea and distance source for bd_carthage ==
Below after some search, but we have to trust the last data in the bd_carthage
{{{
load(paste(datawd,"/dataEDAcommun/chainage/parcours_source/noeud&troncon_final.RData",sep="")) 
write.table(noeuds_final,file=paste(datawd,"/dataEDAcommun/chainage/parcours_source/noeud&troncon_final.csv",sep=""),sep=";",row.names=FALSE,col.names=TRUE)
}}}
{{{
#!sql
drop table noeud_troncon_final;
create table noeud_troncon_final (
Id_BDCARTHAGE varchar,
situation varchar(10), 
distance_mer numeric,
bassin integer,
noeud_mer integer,
niveau integer,
nb_confluence integer,
nb_defluence integer,
noeud_source integer,
distance_source numeric, 
distance_source_cumulee numeric,
statut varchar(30),
strahler integer,
niveau_defluence numeric,
distance_affluent_defluence numeric,
id_carthage_bis varchar(10),
type_noeud_source varchar(30));

COPY noeud_troncon_final from 'C:/base/noeud_troncon_final1.csv' with CSV header delimiter as ';'  NULL as 'NA'
}}}


{{{
#!sql
CREATE INDEX hylcov_node_ID_SOM
  ON hylcov_node
  USING btree
  (ID_SOM);

CREATE INDEX hylcov_arc_id_som_f
  ON hylcov_arc
  USING btree
  (id_som_f);

CREATE INDEX hylcov_arc_ID_SOM
  ON hylcov_node
  USING btree
  (ID_SOM);

CREATE INDEX noeud_troncon_final_Id_BDCARTHAGE
  ON noeud_troncon_final
  USING btree
  (Id_BDCARTHAGE);

drop view if exists hylcov_arc_dist;
create view hylcov_arc_dist as(
select * from hylcov_arc h left join  noeud_troncon_final n on id_som_f=Id_BDCARTHAGE );;
}}}

to shape
{{{
C:\"Program Files"\PostgreSQL\8.4\bin\pgsql2shp -f "C:\eda\exports_shape\hylcov_arc_dist" -p 5433 -u postgres -P postgres -g the_geom -r -k bd_carthage hylcov_arc_dist

}}}
[[Image(source:data/Docs/trac/bd_carthage/distance_source_hylcov_arc.png,600px)]] 

== 3. Get distances into the stationsp2 table ==

=== 31. Saving table stationsp2 from eda20 to bd_carthage ===

{{{
cd C:\eda\backup
C:\"Program Files"\PostgreSQL\8.4\bin\pg_dump  -U postgres -p 5433 -t bd_map.stationsp2 eda2.0> stationsp2.sql
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d bd_carthage -p 5433 -U postgres -f stationsp2.sql
REM une erreur quand il essaye de mettre le schéma bd_map mais là on a pas besoin d'avoir les données dans un schéma
}}}
at this stage it is usefull to check the geometry column table and correct it eventually
=== 32. project stationsp2 again and get table of all projections within 300 ===

Some stations were missing...
 
 * first create two index to speed up queries [[BR]]
bdcarthage 

{{{
#!sql
CREATE INDEX indexhylcov_arc ON hylcov_arc
  USING GIST ( the_geom GIST_GEOMETRY_OPS );
CREATE INDEX indexStationsp2 ON stationsp2
  USING GIST ( the_geom GIST_GEOMETRY_OPS );
}}}


Joining the hylcovarc and stationsp2, the_geom is that of stationsp2
{{{
#!sql
-- creation de la table bd_map_bd_carthage, je rajoute un the_geom a la requete pour pouvoir voir la table dans Qgis  
DROP TABLE IF EXISTS bd_map_bd_carthage;
CREATE TABLE bd_map_bd_carthage as (
     SELECT distinct on (st_codecsp) st_codecsp, distance_source,distance_mer,strahler,id_trhyd, min(distance) as distance, the_geom FROM (
               SELECT st_codecsp, n.distance_source,n.distance_mer,n.strahler ,a.id_trhyd ,CAST(distance(a.the_geom, s.the_geom) as  decimal(15,1)) as distance,s.the_geom
               FROM stationsp2 As s
               INNER JOIN hylcov_arc a ON ST_DWithin(a.the_geom, s.the_geom,300)
		LEFT JOIN noeud_troncon_final n ON a.id_som_f::numeric = n.id_bdcarthage
               WHERE s.the_geom IS NOT NULL
               ORDER BY st_codecsp) as sub
        GROUP BY st_codecsp, distance, distance_source,distance_mer,strahler, id_trhyd, the_geom
 );
alter table  bd_map_bd_carthage add column id serial;
-- mise à jour de la table geometry_columns
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'public', 'bd_map_bd_carthage', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM bd_map_bd_carthage LIMIT 1;

-- creation d'index, clé primaire, et constraintes qui vont bien
alter table bd_map_bd_carthage add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table bd_map_bd_carthage add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'POINT'::text OR the_geom IS NULL);
alter table bd_map_bd_carthage add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
alter table bd_map_bd_carthage ADD CONSTRAINT pk_id PRIMARY KEY(id);
CREATE INDEX indexbd_map_bd_carthage ON bd_map_bd_carthage
  USING GIST ( the_geom GIST_GEOMETRY_OPS );
}}}

=== 33. Saving back table stationsp2 from bd_carthage to eda ===

{{{
cd C:\eda\backup
C:\"Program Files"\PostgreSQL\8.4\bin\pg_dump  -U postgres -p 5433 -t bd_map_bd_carthage bd_carthage> bd_map_bd_carthage.sql
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -p 5433 -U postgres -f bd_map_bd_carthage.sql
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -p 5433 -U postgres -c "alter table bd_map_bd_carthage drop CONSTRAINT  pk_id;"
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -p 5433 -U postgres -c "alter table bd_map_bd_carthage ADD CONSTRAINT pk_id_bd_map_bd_carthage PRIMARY KEY(id);
C:\"Program Files"\PostgreSQL\8.4\bin\psql -d eda2.0 -p 5433 -U postgres -c "alter table bd_map_bd_carthage set schema bd_map"
}}}

== 4. define criteria helping to do the projection ==

Below this table does not only select minimum distance but all ccm segments which are within a 300 m distance from the electrofishing station. [[BR]]
 
{{{
#!sql
DROP TABLE IF EXISTS bd_map.bdmap_ccm2;
CREATE TABLE bd_map.bdmap_ccm2 as (
               SELECT st_codecsp, r.cum_len_sea, r.distance_source, r.strahler as strahlerccm,gid ,CAST(distance(r.the_geom, s.the_geom) as  decimal(15,1)) as distance,s.the_geom 
               FROM bd_map.stationsp2 As s
               INNER JOIN  ccm21.riversegments r ON ST_DWithin(r.the_geom, s.the_geom,300)
               WHERE s.the_geom IS NOT NULL
               ORDER BY st_codecsp 

);
alter table bd_map.bdmap_ccm2 add column id serial;
}}}

The table from the request joining bd_map.bd_map_bd_carthage with the bd_map.bdmap_ccm2 is exported there :
{{{
#!sql
select * from bd_map.bd_map_bd_carthage b join bd_map.bdmap_ccm2 c on b.st_codecsp=c.st_codecsp;
}}}
source:data/dataEDAccm/export_jointure_bdcarthage_ccm.csv [[BR]]
Some processing with R and some results
We are trying to see if distance_source distance_sea and strahler are consistent between bd_carthage and ccm
{{{
# analyse des relations entre distance source et distance mer pour les projections sur la ccm 
library(lattice)
library(ggplot2)
ccm_bd<-read.csv(paste(datawd,"/dataEDAccm/export_jointure_bdcarthage_ccm.csv",sep=""),sep=";",stringsAsFactors = FALSE)
str(ccm_bd)
ccm_bd$distance_source_ccm=ccm_bd$distance_source.1/1000
ccm_bd$cum_len_sea=ccm_bd$cum_len_sea/1000
ccm_bd$pbstrahler=ccm_bd$strahlerccm>(ccm_bd$strahler+1) # which strahler are larger in ccm than in bd_carthage, allowing diff of one 
ratio_dist_mer=ccm_bd$cum_len_sea/ccm_bd$distance_mer
ccm_bd$pb_ratio_dist_mer<-ratio_dist_mer>2&ccm_bd$cum_len_sea>100 # la distance mer est le douible pour la ccm (peut arriver pour les faibles valeurs)
ratio_sources<-ccm_bd$distance_source_ccm/ccm_bd$distance_source
ccm_bd$pbratiosource<-ratio_sources>3 # la distance source ccm est égale à trois fois la distance source bd_carthage
sum(ccm_bd$pbratiosource,na.rm=TRUE)/nrow(ccm_bd) # 12 % de stations mal projettées
sum(ccm_bd$pb_ratio_dist_mer,na.rm=TRUE)/nrow(ccm_bd) # 0.1 % de stations mal projettées
sum(!ccm_bd$pbratiosource&!ccm_bd$pbstrahler&!ccm_bd$pb_ratio_dist_mer,na.rm=TRUE)/nrow(ccm_bd) # on garde 83 % des données
# hitograms
layout(matrix(c(1,2), 1, 2, byrow = TRUE))
hist(ratio_dist_mer[ratio_dist_mer<2],100,main="ratio dist. sea")
abline(v=2,col="red")
mtext(paste(round(sum(ccm_bd$pb_ratio_dist_mer,na.rm=TRUE)/nrow(ccm_bd),3)*100,"% out"),col="red")
hist(ratio_sources[ratio_sources<10],100,main="ratio dist source")
abline(v=3,col="red")
mtext(paste(round(sum(ccm_bd$pbratiosource,na.rm=TRUE)/nrow(ccm_bd),3)*100,"% out"),col="red")
# plots
g<-ggplot(ccm_bd)
g+geom_point(aes(x=distance_source,y=distance_source_ccm,colour=pbstrahler))
g+geom_point(aes(distance_mer,cum_len_sea,colour=pb_ratio_dist_mer))+geom_abline(slope=1,colour="red") +geom_abline(slope=1,intercept=1,colour="green") 
#g+geom_point(aes(distance_mer,cum_len_sea,colour=pbratiosource))+geom_abline(slope=1,colour="red")
g+geom_jitter(aes(x=strahler,y=strahlerccm,colour=pbstrahler),alpha=0.3)
g1<-g+geom_point(aes(x=distance_source,y=distance_source_ccm,colour=pbratiosource))
ccm_bd1<-ccm_bd[ccm_bd$pbstrahler,]
g1+geom_point(aes(x=distance_source,y=distance_source_ccm),size=2,data=ccm_bd1,pch=1)
}}}
__The following criteria have been used to select stations.__ [[BR]]
''' For distance_sea(ccm) > 100 distance_seaccm/distance_sea_bd_carthage>2 ''' [[BR]]
''' Distance_source(ccm)/distance_source(bd_carthage)>3 '''[[BR]]
''' ccm_bd$strahlerccm>(ccm_bd$strahler+1) '''[[BR]]
Finally 83% lines kept but this include doubles [[BR]]
[[Image(source:data/Docs/trac/BDMAP/percentage_out.jpeg,300px)]] 
[[Image(source:data/Docs/trac/BDMAP/comp_dist_source.jpeg,300px)]] 
[[Image(source:data/Docs/trac/BDMAP/comp_dist_mer.jpeg,300px)]] 
[[Image(source:data/Docs/trac/BDMAP/com_strahler.jpeg,300px)]] 

== 5. query to remove the data from  bd_map.bdmap_ccm2 ==
{{{
#!sql

}}}