wiki:CLC Join

Version 96 (modified by cedric, 15 years ago) (diff)

--

This page describes how to calculate surface areas for the different clc variables for one unit catchment

CookBook Eda

this is ticket #56 (detail of code there and in attachment at the bottom of the page)

First an example for Britany

We chose to work on a reduced scale to build the queries as the clc is just damn too big for us, and raster approches are not yet avalaible in postgis, though we fancy using the functions displayed in the developpement page http://trac.osgeo.org/postgis/wiki/WKTRaster.

building a table of clc only for Britany

After several trial neither ST_Crosses nor ST_Contains did work, the example with ST_Intersects provides a surface that is larger than britany, but ST_Crosses missed the area near the coasts as the geom had to be fully within Britany. This is not really a problem as the clipped function in the next paragraph solves this problem and is not using the Britany layer. Howerver it could be usefull to work on a reduced dataset for the clc as the computing time is very high .... A solution will be to build a single area from the catchments layer for the selected country (using the france.wso1_id) and use it in a similar way as we used Britany here....

-- 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_Intersects(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 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);
CREATE INDEX indexclc00_v2_Bretagne ON clc.clc00_v2_Bretagne
  USING GIST ( the_geom GIST_GEOMETRY_OPS ); 

-- run this only once
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;

source:data/Docs/trac/clc/clc_britany.png

Cutting the surface according to the catchments area

Note this one took a long time to run on my computer, one hour for Britany...
It uses ST_Intersects to reduce the search list and ST_Multi(ST_Intersection()) to cut along the primary catchments borders. As a result of using the catchment database only for Britany, the layers that were outside from Britany geographic range are now discarded
Note also that in postgres some lines apear as void of data even when using CAST(ST_AsText() but that those line provide answer FALSE using either ST_IsValid () or ST_IsEmpty() or IS NULL . We had doubts at the beginning as the file was not displayed in Qgis, but the reason was probably not those apparent voids but probably changing ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) AS the_geom into ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) the_geom . Here we say probably as working with Gis still holds some mysteries.

--------------------------------------
--------------------------------------
--DECOUPAGE DES SURFACES
--------------------------------------
--------------------------------------
DROP TABLE IF EXISTS clc.clipped_bretagne;
CREATE TABLE clc.clipped_bretagne AS
SELECT intersected.clcgid, intersected.gid, code_00,the_geom
FROM (SELECT clc.gid as clcgid, c.gid,code_00, 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; --1h12 min
ALTER TABLE clc.clipped_bretagne ADD column id serial PRIMARY KEY;
alter table clc.clipped_bretagne add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.clipped_bretagne add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.clipped_bretagne add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
CREATE INDEX indexclc00clipped_bretagne ON clc.clipped_bretagne
  USING GIST ( the_geom GIST_GEOMETRY_OPS ); 
  -- les limites sont bien meilleures
-- Here to analyse the structure of data in the created table
SELECT ST_AsText(the_geom) from clc.clipped_bretagne limit 20;

source:data/Docs/trac/clc/clipped_bretagne.png
The final file is effectively cut, the following layer was displayed with the purple layer on top of the brown (catchments)
source:data/Docs/trac/coupure des couvertures anthropisées.png

Grouping the surface by gid and code_00

The next steps is to group the different layers, we use the ST_Multi(ST_Collect(()) function to group the POLYGONS into MULTIPOLYGONS (one per line) and before this we do the reverse by transforming all MUTIPOLYGONS into POLYGONS using ST_Dump() . So finally calculating the area using ST_Area is a very quick and simple step. Note the unicity constraint wich ensures that we have only one code_00 per unit catchment.

--------------------------------------
--------------------------------------
--REGROUPEMENT
--------------------------------------
--------------------------------------

DROP TABLE IF EXISTS clc.clipped_bretagne1;
CREATE TABLE clc.clipped_bretagne1 AS (
SELECT gid,code_00,
           ST_Multi(ST_Collect(f.the_geom)) as the_geom
         FROM (SELECT gid, code_00,(ST_Dump(the_geom)).geom As the_geom
                                FROM
                                 clc.clipped_bretagne
                                ) As f
GROUP BY gid,code_00); -- 5s
ALTER TABLE clc.clipped_bretagne1 add column id serial PRIMARY KEY;
alter table clc.clipped_bretagne1 add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.clipped_bretagne1 add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.clipped_bretagne1 add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
CREATE INDEX indexclc00clipped_bretagne1 ON clc.clipped_bretagne1
  USING GIST ( the_geom GIST_GEOMETRY_OPS );
ALTER TABLE clc.clipped_bretagne1 add constraint c_ck_uk  UNIQUE(gid,code_00); -- contrainte d'unicité OK !


ALTER TABLE clc.clipped_bretagne1 add column area numeric;
UPDATE clc.clipped_bretagne1 set area=ST_Area(the_geom); -- 9s

source:data/Docs/trac/clc/clipped_bretagne1.png

So now we have a dataset looking like the following one

SELECT gid,code_00, id,round(area) as area FROM clc.clipped_Bretagne1 order by gid, code_00 limit 10;
gidcode_00idarea( m²)
22512421175989331309
22512423170932890840
225124242106622825609
22512424364891462960
22512431145951081844
225124313198427439
22513311272951413001
22513321199501273065
2251332311022751347
2251332426079720638

Finally we use this dataset to extract data for a new table which will give one value per catchment id. Note that the gid is the gid of the catchment and the table will have to be joined to the riversegments table. We are not using percentage of cover but surface so that we will be able to calculate for any basin, for instance, the cumulated urban area for the whole catchment upstream from one riversegment. We have chosen to group the Corinne Land Cover variables as following

  • 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

Importantly this request will provide repeated lines with the same content it is this necessary to add a DISTINCT ON clause. I had forgotten to include the sum() and group by clause, for instance urban will correspond to the sum of lines whith code_00 111 112 121 122 123 124 131 132 133.

DROP TABLE IF EXISTS clc.surf_area;
CREATE TABLE clc.surf_area AS (
SELECT DISTINCT ON (init.gid) init.gid,
        artificial_surfaces_11_13,
        artificial_vegetated_14,
         arable_land_21,
         permanent_crops_22,
         pastures_23,
         heterogeneous_agricultural_24,
         forest_31,
         natural_32_33,
         wetlands_4,
         inland_waterbodies_51 ,
         marine_water_52
        -- SELECT * 
         FROM (
        SELECT  gid from clc.clipped_bretagne1  ) as init        
        FULL OUTER JOIN (SELECT gid,sum(area) AS artificial_surfaces_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' 
                        GROUP BY gid) AS artificial_surfaces
                       on (init.gid) =(artificial_surfaces.gid)         
        FULL OUTER JOIN (SELECT gid,sum(area) AS artificial_vegetated_14 FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 2)='14'
                        GROUP BY gid) AS artificial_vegetated
                        on artificial_vegetated.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS arable_land_21 FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 2)='21'
                        GROUP BY gid) AS arable_land
                        on arable_land.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS permanent_crops_22 FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 2)='22'
                        GROUP BY gid) AS permanent_crops
                        on permanent_crops.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS pastures_23 FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 2)='23'
                        GROUP BY gid) AS pastures
                        on pastures.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS heterogeneous_agricultural_24 FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 2)='24'
                        GROUP BY gid) AS heterogeneous_agricultural
                        on heterogeneous_agricultural.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS forest_31 FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 2)='31'
                        GROUP BY gid) AS forest
                        ON forest.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(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'
                        GROUP BY gid) AS nature
                        ON nature.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS wetlands_4  FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 1)='4'
                        GROUP BY gid) AS wetlands
                        on wetlands.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS inland_waterbodies_51 FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 2)='51'
                        GROUP BY gid) AS waterbodies
                        on waterbodies.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS marine_water_52 FROM clc.clipped_bretagne1 WHERE 
                        substring(code_00 from 1 for 2)='52'
                        GROUP BY gid) AS marine_water
                        on marine_water.gid =init.gid); --375 ms
ALTER TABLE clc.surf_area ADD CONSTRAINT c_pk_gid_surf_area PRIMARY KEY (gid);
SELECT * FROM clc.surf_area;

The following table will be integrated directy into the riversegments as the gid here corresponds to the riversegment.

DROP TABLE IF EXISTS clc.surf_area1;
CREATE TABLE clc.surf_area1 AS( 
SELECT 
        r.gid,
        area/1e6 as catchment_area,
        CASE WHEN artificial_surfaces_11_13 IS NOT NULL THEN artificial_surfaces_11_13/1e6
        ELSE 0
        END AS artificial_surfaces_11_13,
        CASE WHEN artificial_vegetated_14 IS NOT NULL THEN artificial_vegetated_14/1e6 
        ELSE 0
        END AS artificial_vegetated_14,
        CASE WHEN arable_land_21 IS NOT NULL THEN arable_land_21/1e6 
        ELSE 0
        END AS arable_land_21,
        CASE WHEN permanent_crops_22 IS NOT NULL THEN permanent_crops_22/1e6 
        ELSE 0
        END AS permanent_crops_22,
        CASE WHEN pastures_23 IS NOT NULL THEN pastures_23/1e6 
        ELSE 0
        END AS pastures_23,
        CASE WHEN heterogeneous_agricultural_24 IS NOT NULL THEN heterogeneous_agricultural_24/1e6
        ELSE 0
        END AS heterogeneous_agricultural_24,
        CASE WHEN forest_31 IS NOT NULL THEN forest_31/1e6 
        ELSE 0
        END AS forest_31,
        CASE WHEN natural_32_33 IS NOT NULL THEN natural_32_33/1e6 
        ELSE 0
        END AS natural_32_33,
        CASE WHEN wetlands_4 IS NOT NULL THEN wetlands_4/1e6 
        ELSE 0
        END AS wetlands_4,
        CASE WHEN inland_waterbodies_51 IS NOT NULL THEN inland_waterbodies_51 /1e6 
        ELSE 0
        END AS inland_waterbodies_51,
        CASE WHEN  marine_water_52 IS NOT NULL THEN marine_water_52/1e6 
        ELSE 0
        END AS marine_water_52,
        c.wso1_id,
        c.the_geom      
FROM clc.surf_area p
JOIN ccm21.catchments c ON c.gid=p.gid
JOIN ccm21.riversegments r on r.wso1_id=c.wso1_id
);
SELECT * FROM clc.surf_area1
DROP TABLE IF EXISTS clc.surf_area_analyse;
CREATE TABLE clc.surf_area_analyse AS( 
SELECT 
        gid,
        wso1_id,
        catchment_area,
        artificial_surfaces_11_13+
         artificial_vegetated_14+
         arable_land_21+
         permanent_crops_22+
         pastures_23+
         heterogeneous_agricultural_24+
         forest_31+
         natural_32_33+
         wetlands_4+
         inland_waterbodies_51 +
         marine_water_52 as sum_clc_area ,
        (artificial_surfaces_11_13+
         artificial_vegetated_14+
         arable_land_21+
         permanent_crops_22+
         pastures_23+
         heterogeneous_agricultural_24+
         forest_31+
         natural_32_33+
         wetlands_4+
         inland_waterbodies_51 +
         marine_water_52)/catchment_area AS pourc_clc,
         the_geom
         FROM clc.surf_area1);
ALTER TABLE clc.surf_area_analyse add CONSTRAINT c_pk_gid_area_analyse PRIMARY KEY (gid);
alter table clc.surf_area_analyse add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.surf_area_analyse add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.surf_area_analyse add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
CREATE INDEX indexclc00area_analyse ON clc.surf_area_analyse
  USING GIST ( the_geom GIST_GEOMETRY_OPS );

Proportion between the surface of the clc basins and the catchments in the ccm, in red 0.45 to 0.9 in purple 1.1 à 2 ... The error is acceptable...

source:data/Docs/trac/clc/surf_area_analyse.png

Same but for France

Some slight change to enhance speed in the first request...

--------------------------------------
--------------------------------------
--SURFACE CUT
--------------------------------------
--------------------------------------
DROP TABLE IF EXISTS clc.clipped_france;
CREATE TABLE clc.clipped_france AS
SELECT intersected.clcgid, intersected.gid, code_00,the_geom
FROM (SELECT clc.gid as clcgid, c.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, c.the_geom)) the_geom
        FROM  clc.clc00_v2_europe clc INNER JOIN 
        (select * from ccm21.catchments c WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='France')) as c
        ON  ST_Intersects (c.the_geom,clc.the_geom)
        
       -- AND substring(code_00 from 1 for 1)='1'
       )  AS intersected; --1h12 min
ALTER TABLE clc.clipped_france ADD column id serial PRIMARY KEY;
alter table clc.clipped_france add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.clipped_france add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.clipped_france add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
CREATE INDEX indexclc00clipped_france ON clc.clipped_france
  USING GIST ( the_geom GIST_GEOMETRY_OPS ); 
--------------------------------------
--------------------------------------
--MERGING
--------------------------------------
--------------------------------------
DROP TABLE IF EXISTS clc.clipped_france1;
CREATE TABLE clc.clipped_france1 AS (
SELECT gid,code_00,
           ST_Multi(ST_Collect(f.the_geom)) as the_geom
         FROM (SELECT gid, code_00,(ST_Dump(the_geom)).geom As the_geom
                                FROM
                                 clc.clipped_france
                                ) As f
GROUP BY gid,code_00); -- 5s
ALTER TABLE clc.clipped_france1 add column id serial PRIMARY KEY;
alter table clc.clipped_france1 add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.clipped_france1 add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.clipped_france1 add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
CREATE INDEX indexclc00clipped_france1 ON clc.clipped_france1
  USING GIST ( the_geom GIST_GEOMETRY_OPS );
ALTER TABLE clc.clipped_france1 add constraint c_ck_uk  UNIQUE(gid,code_00); -- contrainte d'unicité
--------------------------------------
--------------------------------------
--AREA
--------------------------------------
--------------------------------------
ALTER TABLE clc.clipped_france1 add column area numeric;
UPDATE clc.clipped_france1 set area=ST_Area(the_geom); -- 9s
--------------------------------------
--------------------------------------
--AREA PER COLUMN FOR CLC TYPE (agregation)
--------------------------------------
--------------------------------------
SELECT gid,code_00, id,round(area) as area FROM clc.clipped_france1 order by gid, code_00 limit 10;
DROP TABLE IF EXISTS clc.surf_area;
CREATE TABLE clc.surf_area AS (
SELECT DISTINCT ON (init.gid) init.gid,
        artificial_surfaces_11_13,
        artificial_vegetated_14,
         arable_land_21,
         permanent_crops_22,
         pastures_23,
         heterogeneous_agricultural_24,
         forest_31,
         natural_32_33,
         wetlands_4,
         inland_waterbodies_51 ,
         marine_water_52
        -- SELECT * 
         FROM (
        SELECT  gid from clc.clipped_france1    ) as init        
        FULL OUTER JOIN (SELECT gid,sum(area) AS artificial_surfaces_11_13 FROM clc.clipped_france1 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' 
                        GROUP BY gid) AS artificial_surfaces
                       on (init.gid) =(artificial_surfaces.gid)         
        FULL OUTER JOIN (SELECT gid,sum(area) AS artificial_vegetated_14 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='14'
                        GROUP BY gid) AS artificial_vegetated
                        on artificial_vegetated.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS arable_land_21 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='21'
                        GROUP BY gid) AS arable_land
                        on arable_land.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS permanent_crops_22 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='22'
                        GROUP BY gid) AS permanent_crops
                        on permanent_crops.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS pastures_23 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='23'
                        GROUP BY gid) AS pastures
                        on pastures.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS heterogeneous_agricultural_24 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='24'
                        GROUP BY gid) AS heterogeneous_agricultural
                        on heterogeneous_agricultural.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS forest_31 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='31'
                        GROUP BY gid) AS forest
                        ON forest.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS natural_32_33 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='32'
                        OR  substring(code_00 from 1 for 2)='33'
                        GROUP BY gid) AS nature
                        ON nature.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS wetlands_4  FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 1)='4'
                        GROUP BY gid) AS wetlands
                        on wetlands.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS inland_waterbodies_51 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='51'
                        GROUP BY gid) AS waterbodies
                        on waterbodies.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS marine_water_52 FROM clc.clipped_france1 WHERE 
                        substring(code_00 from 1 for 2)='52'
                        GROUP BY gid) AS marine_water
                        on marine_water.gid =init.gid); --375 ms
ALTER TABLE clc.surf_area ADD CONSTRAINT c_pk_gid_surf_area PRIMARY KEY (gid);
SELECT * FROM clc.surf_area;
--------------------------------------
--------------------------------------
--REMOVING ZEROS AND JOINING RIVERSEGMENTS AND CATCHMENTS TABLES
--------------------------------------
--------------------------------------
DROP TABLE IF EXISTS clc.surf_area1;
CREATE TABLE clc.surf_area1 AS( 
SELECT 
        r.gid,
        area/1e6 as catchment_area,
        CASE WHEN artificial_surfaces_11_13 IS NOT NULL THEN artificial_surfaces_11_13/1e6
        ELSE 0
        END AS artificial_surfaces_11_13,
        CASE WHEN artificial_vegetated_14 IS NOT NULL THEN artificial_vegetated_14/1e6 
        ELSE 0
        END AS artificial_vegetated_14,
        CASE WHEN arable_land_21 IS NOT NULL THEN arable_land_21/1e6 
        ELSE 0
        END AS arable_land_21,
        CASE WHEN permanent_crops_22 IS NOT NULL THEN permanent_crops_22/1e6 
        ELSE 0
        END AS permanent_crops_22,
        CASE WHEN pastures_23 IS NOT NULL THEN pastures_23/1e6 
        ELSE 0
        END AS pastures_23,
        CASE WHEN heterogeneous_agricultural_24 IS NOT NULL THEN heterogeneous_agricultural_24/1e6
        ELSE 0
        END AS heterogeneous_agricultural_24,
        CASE WHEN forest_31 IS NOT NULL THEN forest_31/1e6 
        ELSE 0
        END AS forest_31,
        CASE WHEN natural_32_33 IS NOT NULL THEN natural_32_33/1e6 
        ELSE 0
        END AS natural_32_33,
        CASE WHEN wetlands_4 IS NOT NULL THEN wetlands_4/1e6 
        ELSE 0
        END AS wetlands_4,
        CASE WHEN inland_waterbodies_51 IS NOT NULL THEN inland_waterbodies_51 /1e6 
        ELSE 0
        END AS inland_waterbodies_51,
        CASE WHEN  marine_water_52 IS NOT NULL THEN marine_water_52/1e6 
        ELSE 0
        END AS marine_water_52,
        c.wso1_id,
        c.the_geom      
FROM clc.surf_area p
JOIN ccm21.catchments c ON c.gid=p.gid
JOIN ccm21.riversegments r on r.wso1_id=c.wso1_id
);
--------------------------------------
--------------------------------------
--COMPARISON OF SURFACES FROM THE CCM AND CLC
--------------------------------------
--------------------------------------
DROP TABLE IF EXISTS clc.surf_area_analyse;
CREATE TABLE clc.surf_area_analyse AS( 
SELECT 
        gid,
        wso1_id,
        catchment_area,
        artificial_surfaces_11_13+
         artificial_vegetated_14+
         arable_land_21+
         permanent_crops_22+
         pastures_23+
         heterogeneous_agricultural_24+
         forest_31+
         natural_32_33+
         wetlands_4+
         inland_waterbodies_51 +
         marine_water_52 as sum_clc_area ,
        (artificial_surfaces_11_13+
         artificial_vegetated_14+
         arable_land_21+
         permanent_crops_22+
         pastures_23+
         heterogeneous_agricultural_24+
         forest_31+
         natural_32_33+
         wetlands_4+
         inland_waterbodies_51 +
         marine_water_52)/catchment_area AS pourc_clc,
         the_geom
         FROM clc.surf_area1);
ALTER TABLE clc.surf_area_analyse add CONSTRAINT c_pk_gid_area_analyse PRIMARY KEY (gid);
alter table clc.surf_area_analyse add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.surf_area_analyse add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.surf_area_analyse add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
CREATE INDEX indexclc00area_analyse ON clc.surf_area_analyse
  USING GIST ( the_geom GIST_GEOMETRY_OPS );

France

DONE calculate for a whole country

Script to create region by region, it sometimes crashes for unexpected reason, after this script there will be missing basins (those at the frontiers of the region). This script takes a long long time, probably nearly a full day of calculation and sometimes it crashes...

/*
"CORSE";1
"LORRAINE";2
"PICARDIE";3
"MIDI-PYRENEES";4
"POITOU-CHARENTES";5
"ALSACE";6
"ILE-DE-FRANCE";7
"NORD-PAS-DE-CALAIS";8
"HAUTE-NORMANDIE";9
"LIMOUSIN";10
"BOURGOGNE";11
"AUVERGNE";12
"RHONE-ALPES";13
"PROVENCE-ALPES-COTE-D'AZUR";14
"FRANCHE-COMTE";15
"BASSE-NORMANDIE";16
"PAYS-DE-LA-LOIRE";17
"AQUITAINE";18
"CENTRE";19
"CHAMPAGNE-ARDENNE";20
"LANGUEDOC-ROUSSILLON";21
"BRETAGNE";22
*/

BEGIN;
DROP TABLE IF EXISTS clc.clipped;
CREATE TABLE clc.clipped AS
SELECT intersected.clcgid, intersected.gid, code_00,the_geom
FROM (SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
        FROM  clc.clc00_v2_europe clc 
        INNER JOIN (
                SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                        SELECT wso1_id FROM ccm21.catchments c 
                        INNER JOIN (SELECT the_geom
                                FROM france.region where reg_id=1) as sub
                        ON ST_Contains(sub.the_geom,c.the_geom))
        )AS sub1
        ON  ST_Intersects (sub1.the_geom,clc.the_geom)
       )  AS intersected;     
ALTER TABLE clc.clipped ADD column id serial PRIMARY KEY;
alter table clc.clipped add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.clipped add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.clipped add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
CREATE INDEX indexclc00clipped ON clc.clipped
  USING GIST ( the_geom GIST_GEOMETRY_OPS ); 
COMMIT;

CREATE INDEX indexclipped
  ON clc.clipped
  USING btree
  (gid);

BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=2) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom))
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;         

BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=3) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                --where c.gid not in (SELECT gid from clc.clipped) -- ralonge trop (on fera le tri à la fin)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;      -- 50 min

-- MIDI-PYRENNEES
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=4) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                --where c.gid not in (SELECT gid from clc.clipped) -- ralonge trop (on fera le tri à la fin)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;      

-- "POITOU-CHARENTES"
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=5) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;   


-- "ALSACE";6
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=6) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;   

-- "ILE-DE-FRANCE";7
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=7) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  

-- "NORD-PAS-DE-CALAIS";8

BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=8) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  

--"HAUTE-NORMANDIE";9

BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=9) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  

--"LIMOUSIN";10

BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=10) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"BOURGOGNE";11
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=11) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  

--"AUVERGNE";12

BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=12) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  

--"RHONE-ALPES";13
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=13) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  

--"PROVENCE-ALPES-COTE-D'AZUR";14
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=14) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"FRANCHE-COMTE";15
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=15) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"BASSE-NORMANDIE";16
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=16) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"PAYS-DE-LA-LOIRE";17
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=17) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"AQUITAINE";18
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=18) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"CENTRE";19
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=19) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"CHAMPAGNE-ARDENNE";20
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=20) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"LANGUEDOC-ROUSSILLON";21
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=21) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  
--"BRETAGNE";22
BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN (
                        SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                                SELECT wso1_id FROM ccm21.catchments c 
                                INNER JOIN (SELECT the_geom
                                        FROM france.region where reg_id=22) as sub
                                ON ST_Contains(sub.the_geom,c.the_geom)
                                )
                )AS sub1
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  

source:data/Docs/trac/clc/clipped.png

CREATE TABLE clc.remaining_catchment AS(        
        SELECT * FROM ccm21.catchments where gid in(
        SELECT gid FROM ccm21.catchments c WHERE c.wso_id IN (SELECT wso_id FROM france.wso WHERE area='France')
        except (select gid from  clc.clipped))
        
);
ALTER TABLE clc.remaining_catchment add CONSTRAINT c_pk_id_remaining_catchment  PRIMARY KEY (gid);
alter table clc.remaining_catchment add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.remaining_catchment add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.remaining_catchment add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);

Using ST_intersection has suppressed the polygons outside catchment

source:data/Docs/trac/clc/clipped_and_catchment.png

Remaining catchmentsis the table of the catchment not already processed

source:data/Docs/trac/clc/remaining_catchment.png

It is a little bit too large to process as one unit

source:data/Docs/trac/clc/remaining_catchment2.png

-- remaining catchments intersecting with france in a table for Qgis
CREATE TABLE clc.remaining_catchment_france AS( 
        SELECT * FROM clc.remaining_catchment  WHERE wso1_id IN (
                SELECT wso1_id FROM clc.remaining_catchment crc JOIN france.departement fd
                ON ST_Intersects(crc.the_geom,fd.the_geom))     
);
ALTER TABLE clc.remaining_catchment_france add CONSTRAINT c_pk_id_remaining_catchment_france  PRIMARY KEY (gid);
alter table clc.remaining_catchment_france add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.remaining_catchment_france add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.remaining_catchment_france add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);

source:data/Docs/trac/clc/remaining_catchment_france.png

Calculating for France and outside from France

BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN 
                (SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                        SELECT wso1_id FROM clc.remaining_catchment_france)  ) AS sub1
                        
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT;  

BEGIN;
INSERT INTO clc.clipped 
        SELECT clc.gid as clcgid, sub1.gid,code_00, ST_Multi(ST_Intersection(clc.the_geom, sub1.the_geom)) the_geom
                FROM  clc.clc00_v2_europe clc 
                INNER JOIN 
                (SELECT gid, c.the_geom FROM ccm21.catchments c where wso1_id IN (
                        SELECT wso1_id FROM clc.remaining_catchment
                        EXCEPT (SELECT wso1_id FROM clc.remaining_catchment_france ))) AS sub1
                        
                ON  ST_Intersects (sub1.the_geom,clc.the_geom);
COMMIT; --~6h45 un plantage     

OOOoooooaaaaaaahhhhhhhhh !!!!!

source:data/Docs/trac/clc/clipped_total.png

--renaming variables for Britany only
ALTER TABLE clc.surf_area RENAME TO surf_area_bretagne;
ALTER TABLE clc.surf_area1 RENAME TO surf_area1_bretagne;
-- grouping rows
DROP TABLE IF EXISTS clc.clipped1;
CREATE TABLE clc.clipped1 AS (
SELECT gid,code_00,
           ST_Multi(ST_Collect(f.the_geom)) as the_geom
         FROM (SELECT gid, code_00,(ST_Dump(the_geom)).geom As the_geom
                                FROM
                                 clc.clipped
                                ) As f
GROUP BY gid,code_00); -- 5s
ALTER TABLE clc.clipped1 add column id serial PRIMARY KEY;
alter table clc.clipped1 add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2);
alter table clc.clipped1 add  CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL);
alter table clc.clipped1 add  CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035);
CREATE INDEX indexclc00clipped1 ON clc.clipped1
  USING GIST ( the_geom GIST_GEOMETRY_OPS );
ALTER TABLE clc.clipped1 add constraint c_ck_uk  UNIQUE(gid,code_00); -- contrainte d'unicité OK !


ALTER TABLE clc.clipped1 add column area numeric;
UPDATE clc.clipped1 set area=ST_Area(the_geom); -- 9s
-- Calculating and grouping surfaces
DROP TABLE IF EXISTS clc.surf_area;
CREATE TABLE clc.surf_area AS (
SELECT DISTINCT ON (init.gid) init.gid,
        artificial_surfaces_11_13,
        artificial_vegetated_14,
         arable_land_21,
         permanent_crops_22,
         pastures_23,
         heterogeneous_agricultural_24,
         forest_31,
         natural_32_33,
         wetlands_4,
         inland_waterbodies_51 ,
         marine_water_52
        -- SELECT * 
         FROM (
        SELECT  gid from clc.clipped1  ) as init        
        FULL OUTER JOIN (SELECT gid,sum(area) AS artificial_surfaces_11_13 FROM clc.clipped1 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' 
                        GROUP BY gid) AS artificial_surfaces
                       on (init.gid) =(artificial_surfaces.gid)         
        FULL OUTER JOIN (SELECT gid,sum(area) AS artificial_vegetated_14 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='14'
                        GROUP BY gid) AS artificial_vegetated
                        on artificial_vegetated.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS arable_land_21 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='21'
                        GROUP BY gid) AS arable_land
                        on arable_land.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS permanent_crops_22 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='22'
                        GROUP BY gid) AS permanent_crops
                        on permanent_crops.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS pastures_23 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='23'
                        GROUP BY gid) AS pastures
                        on pastures.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS heterogeneous_agricultural_24 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='24'
                        GROUP BY gid) AS heterogeneous_agricultural
                        on heterogeneous_agricultural.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS forest_31 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='31'
                        GROUP BY gid) AS forest
                        ON forest.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS natural_32_33 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='32'
                        OR  substring(code_00 from 1 for 2)='33'
                        GROUP BY gid) AS nature
                        ON nature.gid =init.gid
        FULL OUTER JOIN (SELECT gid, sum(area) AS wetlands_4  FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 1)='4'
                        GROUP BY gid) AS wetlands
                        on wetlands.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS inland_waterbodies_51 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='51'
                        GROUP BY gid) AS waterbodies
                        on waterbodies.gid =init.gid
        FULL OUTER JOIN (SELECT gid,sum(area) AS marine_water_52 FROM clc.clipped1 WHERE 
                        substring(code_00 from 1 for 2)='52'
                        GROUP BY gid) AS marine_water
                        on marine_water.gid =init.gid); --375 ms
ALTER TABLE clc.surf_area ADD CONSTRAINT c_pk_gid_surf_area PRIMARY KEY (gid);
DROP TABLE IF EXISTS clc.surf_area_final;
CREATE TABLE clc.surf_area_final AS( 
SELECT 
        r.gid,
        area/1e6 as catchment_area,
        CASE WHEN artificial_surfaces_11_13 IS NOT NULL THEN artificial_surfaces_11_13/1e6
        ELSE 0
        END AS artificial_surfaces_11_13,
        CASE WHEN artificial_vegetated_14 IS NOT NULL THEN artificial_vegetated_14/1e6 
        ELSE 0
        END AS artificial_vegetated_14,
        CASE WHEN arable_land_21 IS NOT NULL THEN arable_land_21/1e6 
        ELSE 0
        END AS arable_land_21,
        CASE WHEN permanent_crops_22 IS NOT NULL THEN permanent_crops_22/1e6 
        ELSE 0
        END AS permanent_crops_22,
        CASE WHEN pastures_23 IS NOT NULL THEN pastures_23/1e6 
        ELSE 0
        END AS pastures_23,
        CASE WHEN heterogeneous_agricultural_24 IS NOT NULL THEN heterogeneous_agricultural_24/1e6
        ELSE 0
        END AS heterogeneous_agricultural_24,
        CASE WHEN forest_31 IS NOT NULL THEN forest_31/1e6 
        ELSE 0
        END AS forest_31,
        CASE WHEN natural_32_33 IS NOT NULL THEN natural_32_33/1e6 
        ELSE 0
        END AS natural_32_33,
        CASE WHEN wetlands_4 IS NOT NULL THEN wetlands_4/1e6 
        ELSE 0
        END AS wetlands_4,
        CASE WHEN inland_waterbodies_51 IS NOT NULL THEN inland_waterbodies_51 /1e6 
        ELSE 0
        END AS inland_waterbodies_51,
        CASE WHEN  marine_water_52 IS NOT NULL THEN marine_water_52/1e6 
        ELSE 0
        END AS marine_water_52,
        c.wso1_id,
        c.the_geom      
FROM clc.surf_area p
JOIN ccm21.catchments c ON c.gid=p.gid
JOIN ccm21.riversegments r on r.wso1_id=c.wso1_id
);

saving the result of the request

C:\"Program Files"\PostgreSQL\8.4\bin\pg_dump.exe --host localhost --port 5433 --username postgres --format custom --blobs --verbose --table="clc.clipped" --file "C:\Documents and Settings\cedric\Mes documents\Migrateur\eda\CLC\clipped.backup" "eda2.0"
C:\"Program Files"\PostgreSQL\8.4\bin\pg_dump.exe --host localhost --port 5433 --username postgres --format custom --blobs --verbose --table="clc.clipped1" --file "C:\Documents and Settings\cedric\Mes documents\Migrateur\eda\CLC\clipped1.backup" "eda2.0"
C:\"Program Files"\PostgreSQL\8.4\bin\pg_dump.exe --host localhost --port 5433 --username postgres --format custom --blobs --verbose --table="clc.surf_area" --file "C:\Documents and Settings\cedric\Mes documents\Migrateur\eda\CLC\surf_area.backup" "eda2.0"

TODO integrate the data in R

Attachments (3)

Download all attachments as: .zip