before i run this code, we need to create a clc schema within eda2 on my postgres. right click on schema and choose new schema
then use the following code in the cmd window to load the clc data from the c drive to the database
c: ---to change from h drive to c drive cd data ---to locate thr data folder psql -U postgres -f "clc00_v2_europe.sql" eda2 ---to get the data from the file eda2.clc into postgres in the eda2 database --open pgadmin and run this sql qyery to check its right select count(*) from clc.clc00_v2_europe --20709 good --ctrl c will cancel a command if you know you have messed it up --if something goes into the wrong schema - for example if you havent set up the clc schema, or the.sql file doenst -- --know where to out it, you can use this command to rest its location
--alter table mytable set schema myschema --e.g. alter table clc00_v2_europe set schema clc
riverbasin districts from EU, download them at my ftp/cedric/elvira OPEN COMMAND batch cmd CREATE SCHEMA psql -U postgres -c "create schema european_wise2008" eda2 RESTORE TABLE psql -U postgres -f "european_wise2008.rbd_f1v3.sql" eda2 select count(*) from european_wise2008.rbd_f1v3--215
-- Irish wrbd select name_engl, the_geom from european_wise2008.rbd_f1v3 where gid in (30,31,32,33,34,35,36,37,38,39); -- 39 is already inserted -- what is the sructure of the output table ? --"Neagh Bann" --"North Western" --"South Eastern" --"Shannon" --"South Western" --"Western" --"Neagh Bann" --"North Western" --"North Eastern" --"Eastern" select * from europe.wso limit 10;-- id, wso_id, area -- I just need the two second columns, below a request spatial joining riversegments and the wise layer select distinct on (wso_id) wso_id, name_engl as area from ccm21.riversegments r join (SELECT the_geom, name_engl FROM european_wise2008.rbd_f1v3 As f where gid in (30,31,32,33,34,35,36,37,38) ) as sub ON ST_Intersects(sub.the_geom,r.the_geom) limit 10; --83746;"Neagh Bann" --83747;"Shannon" --83749;"South Eastern" --83772;"North Western" --83773;"North Western" --83783;"Eastern" --83787;"Eastern" --83798;"South Western" --83918;"North Western" --83932;"North Western" -- i'll insert this in europe.wso insert into europe.wso(wso_id,area) select distinct on (wso_id) wso_id, name_engl as area from ccm21.riversegments r join (SELECT the_geom, name_engl FROM european_wise2008.rbd_f1v3 As f where gid in (30,31,32,33,34,35,36,37,38) ) as sub ON ST_Intersects(sub.the_geom,r.the_geom) ;--472 -- Which basin names now ? select max(id) from europe.wso; --8419 select distinct on (area) area from (select * from europe.wso where id >(8419-472))sub; /* "Eastern" "Neagh Bann" "North Eastern" "North Western" "Shannon" "South Eastern" "South Western" */
Dumping and reloading the files already done for western RBD
I'm not sure you'll need it. I suggest doing the work for all basins in Ireland including western, anyway here it is
-- saving for Elvira, the file is named clc_western_tables.zip on my ftp pg_dump -U postgres -h 192.168.1.104 -f "clc_western_tables.sql" --table clc.clipped_western --table clc.surf_area_western --table clc.surf_area_western_final --verbose eda2 -- reloading psql -U postgres -f "clc_western_tables.sql" eda2
Checking data
select count(*) from clc.surf_area_western_final --958 yes ok
-------------------------------------------- --------------------------------------------- -- Corinne Landcover --------------------------------------------- --------------------------------------------- CREATE INDEX indexcatchment_sp ON ccm21.catchments USING GIST ( the_geom GIST_GEOMETRY_OPS ); -------------------------------------- /* -- Elvira : here I'm renaming the tables labelled something_ireland to something_western -- so that you can run the code below alter table clc.clipped_ireland rename to clipped_western; alter table clc.clipped_ireland1 rename to clipped_western1; alter table clc.surf_area_ireland rename to surf_area_western; alter table clc.surf_area_ireland_final rename to surf_area_western_final; */ /* --have switched these files around to relabel something_western to something_irleand alter table clc.clipped_western rename to clipped_ireland;--renamed fine alter table clc.clipped_western1 rename to clipped_ireland1;---this one doesnt exist - do i need it? alter table clc.surf_area_western rename to surf_area_ireland;---renamed fine alter table clc.surf_area_western_final rename to surf_area_ireland_final;---this one doesnt exist - although it looks like it is there - but the number of rows is 0i did do the count above and it worked though */ -------- --SURFACE CUT -------------------------------------- -------------------------------------- -- Below i'm using a BEGIN COMMIT statement so that if the program drops at some point it will not loose everything -- Eastern 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 WHERE c.wso_id IN (SELECT wso_id FROM europe.wso WHERE area='Eastern') ) )AS sub1 ON ST_Intersects (sub1.the_geom,clc.the_geom) ) AS intersected; COMMIT; -- launch first to there
BEGIN; --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); ---all fine to here /* "Eastern"---done 5059 rows ----Query returned successfully with no result in 33812360 ms. "Western"---- done 14709 rows retrieved.Query returned successfully with no result in 23451828 ms. "Neagh Bann"----done 19007 rows. Query returned successfully with no result in 47476391 ms. "North Eastern"--- one 20542 rows. Query returned successfully with no result in 8935422 ms. "North Western"----done 30667 rows. Query returned successfully with no result in 48430063 ms. "Shannon"----done 47274 rows. Query returned successfully with no result in 102877984 ms. "South Eastern"----57130 rows. Query returned successfully with no result in 51168469 ms. "South Western" ---66480 rows all done. Query returned successfully with no result in 37618156 ms. */ -- Neagh Bann 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 WHERE c.wso_id IN (SELECT wso_id FROM europe.wso WHERE area='Neagh Bann') ) )AS sub1 ON ST_Intersects (sub1.the_geom,clc.the_geom); COMMIT; -- add the others
This all seemes to work OK,hve clipped the clc for all the rbds (8 in total)
-------------------------------------- -------------------------------------- --SOME CATCHMENTS ARE MISSING -------------------------------------- -------------------------------------- DROP TABLE if exists clc.remaining_catchments; CREATE TABLE clc.remaining_catchments AS( SELECT gid,wso1_id,wso_id,the_geom FROM ccm21.catchments where gid in( select c.gid from wrbd.outline_of_ireland oi join ccm21.catchments c on st_intersects(oi.the_geom,c.the_geom) except select gid from clc.clipped )); ALTER TABLE clc.remaining_catchments add CONSTRAINT c_pk_id_remaining_catchment PRIMARY KEY (gid); alter table clc.remaining_catchments add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2); alter table clc.remaining_catchments add CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL); alter table clc.remaining_catchments add CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035); -- Now we will try to insert the latest CREATE INDEX indexremainingcatchments ON clc.remaining_catchments USING GIST ( the_geom GIST_GEOMETRY_OPS ); 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_catchments) )AS sub1 ON ST_Intersects (sub1.the_geom,clc.the_geom); -------------------------------------- -------------------------------------- --MERGING -------------------------------------- -------------------------------------- -- at this stage I would like to see the result in Qgis, are there holes ? -- from now on I'm using the name "ireland" to indicate that everything is Ok DROP TABLE IF EXISTS clc.clipped_ireland1; CREATE TABLE clc.clipped_ireland1 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);--done - 7th march EdeE. looks ok - desnt seem to be any holes ALTER TABLE clc.clipped_ireland1 add column id serial PRIMARY KEY; alter table clc.clipped_ireland1 add CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2); alter table clc.clipped_ireland1 add CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL); alter table clc.clipped_ireland1 add CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3035); CREATE INDEX indexclc00clipped_ireland1 ON clc.clipped_ireland1 USING GIST ( the_geom GIST_GEOMETRY_OPS ); ALTER TABLE clc.clipped_ireland1 add constraint c_ck_uk_ireland UNIQUE(gid,code_00); -- contrainte d'unicité --------------------------------------done 7th march edeE -------------------------------------- --AREA -------------------------------------- -------------------------------------- -- what I'm doing below is transforming from a table having one line per clc type and per unit basin (where I have the type id and a geometry) into one line, where I calculate the total surface of each landcover category. ALTER TABLE clc.clipped_ireland1 add column area numeric; UPDATE clc.clipped_ireland1 set area=ST_Area(the_geom); ------------------------done 7th march ede E-------------- -------------------------------------- --AREA PER COLUMN FOR CLC TYPE (agregation) -------------------------------------- -------------------------------------- SELECT gid,code_00, id,round(area) as area FROM clc.clipped_ireland1 order by gid, code_00 limit 10; DROP TABLE IF EXISTS clc.surf_area_ireland; CREATE TABLE clc.surf_area_ireland 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_ireland1 ) as init FULL OUTER JOIN (SELECT gid,sum(area) AS artificial_surfaces_11_13 FROM clc.clipped_ireland1 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_ireland1 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_ireland1 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_ireland1 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_ireland1 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_ireland1 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_ireland1 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_ireland1 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_ireland1 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_ireland1 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_ireland1 WHERE substring(code_00 from 1 for 2)='52' GROUP BY gid) AS marine_water on marine_water.gid =init.gid); --3406ms EdeE ALTER TABLE clc.surf_area_ireland ADD CONSTRAINT c_pk_gid_surf_area_ireland PRIMARY KEY (gid); SELECT * FROM clc.surf_area_ireland; -------------------------------------- -------------------------------------- --REMOVING ZEROS AND JOINING RIVERSEGMENTS AND CATCHMENTS TABLES -------------------------------------- -------------------------------------- -- this table drops the previous one but final calculations are stored in surf_area_final DROP TABLE IF EXISTS clc.surf_area_ireland_final; CREATE TABLE clc.surf_area_ireland_final AS( SELECT c.gid, C.area/1e6 as catchment_area, CASE WHEN p.artificial_surfaces_11_13 IS NOT NULL THEN p.artificial_surfaces_11_13/1e6 ELSE 0 END AS artificial_surfaces_11_13, CASE WHEN p.artificial_vegetated_14 IS NOT NULL THEN p.artificial_vegetated_14/1e6 ELSE 0 END AS artificial_vegetated_14, CASE WHEN p.arable_land_21 IS NOT NULL THEN p.arable_land_21/1e6 ELSE 0 END AS arable_land_21, CASE WHEN p.permanent_crops_22 IS NOT NULL THEN p.permanent_crops_22/1e6 ELSE 0 END AS permanent_crops_22, CASE WHEN p.pastures_23 IS NOT NULL THEN p.pastures_23/1e6 ELSE 0 END AS pastures_23, CASE WHEN p.heterogeneous_agricultural_24 IS NOT NULL THEN p.heterogeneous_agricultural_24/1e6 ELSE 0 END AS heterogeneous_agricultural_24, CASE WHEN p.forest_31 IS NOT NULL THEN p.forest_31/1e6 ELSE 0 END AS forest_31, CASE WHEN p.natural_32_33 IS NOT NULL THEN p.natural_32_33/1e6 ELSE 0 END AS natural_32_33, CASE WHEN p.wetlands_4 IS NOT NULL THEN p.wetlands_4/1e6 ELSE 0 END AS wetlands_4, CASE WHEN p.inland_waterbodies_51 IS NOT NULL THEN p.inland_waterbodies_51 /1e6 ELSE 0 END AS inland_waterbodies_51, CASE WHEN p.marine_water_52 IS NOT NULL THEN p.marine_water_52/1e6 ELSE 0 END AS marine_water_52, c.wso1_id, c.the_geom FROM clc.surf_area_ireland p JOIN ccm21.catchments c ON c.gid=p.gid JOIN ccm21.riversegments r on r.wso1_id=c.wso1_id );
This is fine now although cedric changed the code above a little bit and i forgot to save it Seems that the total areas caluclated from the clc database are larger in a lot of cases than the surface area of the ccm catchments (is this something to do with the lakes/)
working on some code to express the areas of land use as percentages of the sum of all clc land uses (rather than ccm surface area)
--this bit of code sums up all the clc areas and comapres them to the ccm surface areas DROP TABLE IF EXISTS clc.surf_area_analyse; CREATE TABLE clc.surf_area_analyse AS( SELECT 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_area_ireland_final);
--this bit of code divides each categrpry by the sum of all clc categories to get percent land cover. does the clc categorgy wetlands mean bog? i think so DROP TABLE IF EXISTS clc.percent; CREATE TABLE clc.percent AS( SELECT 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_of_CLC_areas, artificial_surfaces_11_13/(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 pourc_artifical_surfaces, artificial_vegetated_14/(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 pourc_artificial_vegetated, arable_land_21/(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 pourc_arable_land, permanent_crops_22/(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 pourc_permanent_crops, pastures_23/(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 pourc_pastures, heterogeneous_agricultural_24/(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 pourc_heterogeneous_agricultural, forest_31/(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 pourc_forest, natural_32_33/(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 pourc_natural, wetlands_4/(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 pourc_wetlands, inland_waterbodies_51/(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 pourc_inland_waterbodies, marine_water_52/(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 pourc_marine_water, the_geom FROM clc.surf_area_ireland_final);
---this bit of code lets you takestuff (clc.surf_area_ireland_final) out of pgadmin and save it to my c drive C:\data>pg_dump -U postgres -f "clc_final.sql" --table clc.surf_area_ireland_fin al --verbose eda2
The gid's are needed to work the ccm21.upstream(_gid) function. Here is how we integrate them in the table
alter table clc.surf_area_ireland_final add column gid integer; alter table clc.surf_area_ireland_final add constraint c_uk_gid unique(gid); update clc.surf_area_ireland_final set gid=riversegments.gid from ccm21.riversegments where riversegments.wso1_id=surf_area_ireland_final.wso1_id; --9064 select count(*) from clc.surf_area_ireland_final; --10035
Last modified 13 years ago
Last modified on May 17, 2012 9:04:36 PM
