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 attached file)

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;

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

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

The final file is effectively cut, the following layer was displayed with the purple layer on top of the brown (catchments)

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

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;

gid	code_00	id	area( m²)
225124	211	7598	9331309
225124	231	7093	2890840
225124	242	10662	2825609
225124	243	6489	1462960
225124	311	4595	1081844
225124	313	198	427439
225133	112	7295	1413001
225133	211	9950	1273065
225133	231	1022	751347
225133	242	6079	720638

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

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

A whole country

TODO calculate for a whole country

TODO integrate the data in R