Integrating Corine Land Cover data into Postgre

[celine]

-- Un essai sur la Bretagne pour voir
DROP TABLE IF EXISTS clc.clipped_bretagne;
CREATE TABLE clc.clipped_bretagne AS
SELECT intersected.clcgid, intersected.gid, clipped_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) AS clipped_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1')  AS intersected;--3721 lines   
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;
-- ajout dans la table geometry_columns pour référencement rapide sous Qgis
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'clc', 'clipped_bretagne', 'the_geom', ST_CoordDim(clipped_geom), ST_SRID(clipped_geom), MULTIPOLYGON
FROM clc.clipped_bretagne  LIMIT 1;

-- il ne s'affiche pas... essai réduit OK
DROP TABLE IF EXISTS clc.clipped_bretagne;
CREATE TABLE clc.clipped_bretagne AS
SELECT clc.gid as clcgid,c.gid, clc.the_geom FROM clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1';
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;

-- la requete globale
SELECT intersected.clcgid, intersected.gid, clipped_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Intersection(clc.the_geom, c.the_geom) AS clipped_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        where substring(code_00 from 1 for 1)='1' )  AS intersected;            

-- script de départ (modifié ensuite)
CREATE TABLE clc.clipped_bretagne AS
SELECT intersected.clcgid, intersected.gid, clipped_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) clipped_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1')  AS intersected;--3721 lines  
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;

-3721 lines  
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;
-- ajout dans la table geometry_columns pour référencement rapide sous Qgis
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'clc', 'clipped_bretagne', 'the_geom', ST_CoordDim(clipped_geom), ST_SRID(clipped_geom), MULTIPOLYGON
FROM clc.clipped_bretagne  LIMIT 1;

-- ci dessous j'essaye d'analyser la structure interne des données pour comprendre quels sont les problèmes

SELECT intersected.clcgid, intersected.gid, intersected.the_geom
FROM (SELECT clc.gid as clcgid, c.gid, CAST(ST_AsText(ST_Multi(ST_Intersection(clc.the_geom, c.the_geom))) as VARCHAR(150)) as the_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1')  AS intersected;


-- Pour comparaison (ressemble beaucoup)

SELECT CAST(ST_AsText(the_geom) as varchar(150)) from clc.clc00_v2_europe clc limit 20;


-- Je rééssaye
DROP TABLE IF EXISTS clc.clipped_bretagne;
CREATE TABLE clc.clipped_bretagne AS
SELECT intersected.clcgid, intersected.gid, the_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) the_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1')  AS intersected;--3721 lines  
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;
-- Here to analyse the structure of data in the created table
SELECT ST_AsText(the_geom) from clc.clipped_bretagne limit 20;

-- celui là marche sur Postgis (et je ne comprends pas ce qui ne marchait pas peut être le as the_geom a la place de the_geom, ou alors mon Qgis....)

-- I'm adding some original data tables restricted for an area for Qgis easy use
-- extraction of a clc table for Britany
DROP TABLE IF EXISTS clc.clc00_v2_Bretagne;
CREATE TABLE clc.clc00_v2_Bretagne AS
SELECT * FROM clc.clc00_v2_europe where gid IN (
SELECT gid FROM  clc.clc00_v2_europe clc JOIN
    (SELECT the_geom FROM france.region where code_reg='53') as sub
    ON ST_Contains(sub.the_geom,clc.the_geom));
ALTER TABLE clc.clc00_v2_Bretagne ADD CONSTRAINT c_pk_gid PRIMARY KEY  (gid);
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'clc', 'clc00_v2_Bretagne', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM clc.clc00_v2_europe LIMIT 1;

-- Creating a table for riversegments Bretagne

DROP TABLE IF EXISTS ccm21.riversegments_Bretagne;   
CREATE TABLE  ccm21.riversegments_Bretagne AS
    SELECT * FROM ccm21.riversegments
        WHERE wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne');
ALTER TABLE ccm21.riversegments_Bretagne ADD CONSTRAINT c_pk_gid_riversegments_Bretagne PRIMARY KEY  (gid);
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'ccm21', 'riversegments_Bretagne', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM ccm21.riversegments_Bretagne LIMIT 1;

-- creating a table for catchments Bretagne.
DROP TABLE IF EXISTS ccm21.catchments_Bretagne;   
CREATE TABLE  ccm21.catchments_Bretagne AS
    SELECT * FROM ccm21.catchments
        WHERE wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne');  
ALTER TABLE ccm21.catchments_Bretagne ADD CONSTRAINT c_pk_gid_catchments_Bretagne PRIMARY KEY  (gid);
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'ccm21', 'catchments_Bretagne', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM ccm21.catchments_Bretagne LIMIT 1;

[cedric]

ST_Area - Returns the area of the surface if it is a polygon or multi-polygon. For "geometry" type area is in SRID units. For "geography" area is in square meters.

3035
+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

[cedric]

-- Un essai sur la Bretagne pour voir
DROP TABLE IF EXISTS clc.clipped_bretagne;
CREATE TABLE clc.clipped_bretagne AS
SELECT intersected.clcgid, intersected.gid, clipped_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) AS clipped_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1')  AS intersected;--3721 lines   
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;
-- ajout dans la table geometry_columns pour référencement rapide sous Qgis
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'clc', 'clipped_bretagne', 'the_geom', ST_CoordDim(clipped_geom), ST_SRID(clipped_geom), MULTIPOLYGON
FROM clc.clipped_bretagne  LIMIT 1;

-- il ne s'affiche pas... essai réduit OK
DROP TABLE IF EXISTS clc.clipped_bretagne;
CREATE TABLE clc.clipped_bretagne AS
SELECT clc.gid as clcgid,c.gid, clc.the_geom FROM clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1';
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;

-- la requete globale
SELECT intersected.clcgid, intersected.gid, clipped_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Intersection(clc.the_geom, c.the_geom) AS clipped_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        where substring(code_00 from 1 for 1)='1' )  AS intersected;            

-- script de départ (modifié ensuite)
CREATE TABLE clc.clipped_bretagne AS
SELECT intersected.clcgid, intersected.gid, clipped_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) clipped_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1')  AS intersected;--3721 lines  
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;

-3721 lines  
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;
-- ajout dans la table geometry_columns pour référencement rapide sous Qgis
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'clc', 'clipped_bretagne', 'the_geom', ST_CoordDim(clipped_geom), ST_SRID(clipped_geom), MULTIPOLYGON
FROM clc.clipped_bretagne  LIMIT 1;

-- ci dessous j'essaye d'analyser la structure interne des données pour comprendre quels sont les problèmes

SELECT intersected.clcgid, intersected.gid, intersected.the_geom
FROM (SELECT clc.gid as clcgid, c.gid, CAST(ST_AsText(ST_Multi(ST_Intersection(clc.the_geom, c.the_geom))) as VARCHAR(150)) as the_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1')  AS intersected;


-- Pour comparaison (ressemble beaucoup)

SELECT CAST(ST_AsText(the_geom) as varchar(150)) from clc.clc00_v2_europe clc limit 20;


-- Je rééssaye
DROP TABLE IF EXISTS clc.clipped_bretagne;
CREATE TABLE clc.clipped_bretagne AS
SELECT intersected.clcgid, intersected.gid, the_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) the_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN ccm21.catchments c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne')
        AND substring(code_00 from 1 for 1)='1')  AS intersected;--3721 lines  
ALTER TABLE clc.clipped_bretagne add column id serial PRIMARY KEY;
-- Here to analyse the structure of data in the created table
SELECT ST_AsText(the_geom) from clc.clipped_bretagne limit 20;

-- celui là marche sur Postgis (et je ne comprends pas ce qui ne marchait pas peut être le as the_geom a la place de the_geom, ou alors mon Qgis....)

-- I'm adding some original data tables restricted for an area for Qgis easy use
-- extraction of a clc table for Britany
DROP TABLE IF EXISTS clc.clc00_v2_Bretagne;
CREATE TABLE clc.clc00_v2_Bretagne AS
SELECT * FROM clc.clc00_v2_europe where gid IN (
SELECT gid FROM  clc.clc00_v2_europe clc JOIN
    (SELECT the_geom FROM france.region where code_reg='53') as sub
    ON ST_Crosses(sub.the_geom,clc.the_geom));
    
ALTER TABLE clc.clc00_v2_Bretagne ADD CONSTRAINT c_pk_gid PRIMARY KEY  (gid);
alter table clc.clc00_v2_Bretagne add column id serial PRIMARY KEY;
alter table clc.clc00_v2_Bretagne add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.clc00_v2_Bretagne add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.clc00_v2_Bretagne add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);


INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'clc', 'clc00_v2_Bretagne', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM clc.clc00_v2_europe LIMIT 1;


CREATE INDEX indexclc00_v2_Bretagne ON clc.clc00_v2_Bretagne
  USING GIST ( the_geom GIST_GEOMETRY_OPS ); 


-- Creating a table for riversegments Bretagne

DROP TABLE IF EXISTS ccm21.riversegments_Bretagne;   
CREATE TABLE  ccm21.riversegments_Bretagne AS
    SELECT * FROM ccm21.riversegments
        WHERE wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne');
ALTER TABLE ccm21.riversegments_Bretagne ADD CONSTRAINT c_pk_gid_riversegments_Bretagne PRIMARY KEY  (gid);
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'ccm21', 'riversegments_Bretagne', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM ccm21.riversegments_Bretagne LIMIT 1;

-- creating a table for catchments Bretagne.
DROP TABLE IF EXISTS ccm21.catchments_Bretagne;   
CREATE TABLE  ccm21.catchments_Bretagne AS
    SELECT * FROM ccm21.catchments
        WHERE wso_id IN (SELECT wso_id FROM france.wso WHERE area='Bretagne');  
ALTER TABLE ccm21.catchments_Bretagne ADD CONSTRAINT c_pk_gid_catchments_Bretagne PRIMARY KEY  (gid);
INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'ccm21', 'catchments_Bretagne', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), GeometryType(the_geom)
FROM ccm21.catchments_Bretagne LIMIT 1;



SELECT * FROM clc.clipped_bretagne;

SELECT SUM(ST_Area(the_geom)) as area, gid, FROM clc.clipped_bretagne;

SELECT ST_AsText(the_geom) from clc.clipped_bretagne limit 20;
--MULTIPOLYGONS
/*
Unfortunately geometry collections are not well-supported by GIS tools. 
To prevent ST_Collect from returning a Geometry Collection when collecting MULTI geometries,
 one can use the below trick that utilizes ST_Dump  to expand the MULTIs out to singles and then regroup them.
 */




 
DROP TABLE IF EXISTS clc.clipped_bretagne1;
CREATE TABLE clc.clipped_bretagne1 AS (
SELECT gid,
	   ST_Multi(ST_Collect(f.the_geom)) as the_geom
	 FROM (SELECT gid, (ST_Dump(the_geom)).geom As the_geom
				FROM
				 clc.clipped_bretagne
				) As f
GROUP BY gid);
ALTER TABLE clc.clipped_bretagne1 add column id serial PRIMARY KEY;

SELECT ST_AsText(the_geom) from clc.clipped_bretagne1 ;

INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'clc', 'clipped_bretagne1', 'the_geom', ST_CoordDim(the_geom), ST_SRID(the_geom), 'MULTIPOLYGON'
FROM clc.clipped_bretagne1  LIMIT 1;

-- j'ai des trous

SELECT ST_AsText(the_geom) FROM (SELECT gid, (ST_Dump(the_geom)).geom As the_geom
				FROM
				 clc.clipped_bretagne) AS f

-- J'ai bien des trous mais the geom is null ne marche pas 4889

SELECT ST_IsValid (the_geom) FROM (SELECT gid, (ST_Dump(the_geom)).geom As the_geom
				FROM
				 clc.clipped_bretagne) AS f

				 -- TRUE partout

SELECT ST_IsEmpty(the_geom) FROM (SELECT gid, (ST_Dump(the_geom)).geom As the_geom
				FROM
				 clc.clipped_bretagne) AS f;
				 -- FALSE partout...; ????

--* bon j'essaye dans Qgis

select count(*) as count, code_00  from clc.clc00_v2_europe group by code_00 order by code_00; -- 

Delete from clc.clc00_v2_europe where code_00='121'

[cedric]

-- here I'm keeping the surface so that we will b able to calculate for any basin
-- or group of basins a percentage of surface according to the basin surface
DROP TABLE IF EXISTS clc.pourc_area;
CREATE TABLE clc.pourc_area AS (
SELECT urban.gid,

urban_11_13,

green_urban_14,
arable_land_21,
plantations_22,
pastures_23,
crops_natural_24,
forest_31,
natural_32_33,
wetlands_4,
water_51 ,
seawater_52
FROM (

SELECT distinct ON (gid) gid from clc.clipped_bretagne1 ) as init
FULL OUTER JOIN (SELECT gid,area AS urban_11_13 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='11'
OR substring(code_00 from 1 for 2)='12'
OR substring(code_00 from 1 for 2)='13' ) AS urban

on init.gid =urban.gid

FULL OUTER JOIN (SELECT gid,area AS green_urban_14 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='14') AS green_urban
on green_urban.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS arable_land_21 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='21') AS arable_land
on arable_land.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS plantations_22 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='22') AS plantations
on plantations.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS pastures_23 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='23') AS pastures
on pastures.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS crops_natural_24 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='24') AS crops_natural
on crops_natural.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS forest_31 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='31') AS forest
on forest.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS natural_32_33 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='32'
OR substring(code_00 from 1 for 2)='33') AS nature
on nature.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS wetlands_4 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 1)='4') AS wetlands
on wetlands.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS water_51 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='51') AS waterbodies
on waterbodies.gid =init.gid

FULL OUTER JOIN (SELECT gid,area AS seawater_52 FROM clc.clipped_bretagne1 WHERE

substring(code_00 from 1 for 2)='52') AS seawater
on seawater.gid =init.gid); --375 ms

[cedric]

See file attachment for draft

[cedric]

Milestone Data integration deleted