CLC Join: ticket#56.sql

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

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;


-- Un essai sur la Bretagne pour voir


-- la requete globale
SELECT intersected.clcgid, intersected.gid, the_geom
FROM (SELECT clc.gid as clcgid, c.gid, ST_Intersection(clc.the_geom, c.the_geom) AS the_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




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
-- ca marche a priori pas de trous
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'

select * from clc.clc00_v2_Bretagne
--------------------------------------
--------------------------------------
--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;
--------------------------------------
--------------------------------------
--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

SELECT gid,code_00, id,round(area) as area FROM clc.clipped_Bretagne1 order by gid, code_00 limit 10;
-- 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.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;


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