1

I have two vector layers (1) with 107 subdistricts (2) soil map with classification of the soil value. I like to calculate the share of the area of the soil map with values classified as "very_high" for each subdistrict. I would like to calculate it with postgresql code. I tried it with

`select g.gid, g.name as subdistrict, round(ST_area(ST_union(ST_INTERSECTION(g.geom,b.geom)))::numeric/100/g.area_ha, 2) as "share_value_very_high",
ST_union(ST_INTERSECTION(g.geom,b.geom)) as geom_union
FROM soil.soilmap AS b, structure.subdistrict as g
WHERE ST_Intersects(g.geom, b.geom) and b.value = 'very_high'
group by g.gid, g.name, b.value;`

It works well and the calcualtion was fine but the resulting table contains only those subdistricts , all in all 78, having soils of very high value. So I tried a left join to get all subdistricts and join the soil map as following

select g.gid, g.name as subdistrict, round(ST_area(ST_union(ST_INTERSECTION(g.geom,b.geom)))::numeric/100/g.area_ha, 2) as "share_value_very_high",
ST_union(ST_INTERSECTION(g.geom,b.geom)) as geom_union
FROM structure.subdistrict as g
left join soil.soilmap AS b 
on ST_within ((b.geom, g.geom) 
WHERE b.value = 'very_high'
group by g.gid, g.name, b.value;

But then the calculation was completely wrong e.g. while the correct share would be 90 % the resulting calculation was only 1 %. Also, I have now only 66 subdistricts. Obviously, I have done a big mistake but where is it?.

  • when I substitute 'ST_within' with 'ST_intersects' in the latter code the calculation of area was correct again but still subdistricts without soil values of very high value doesn´t appear – oppia04 Sep 26 '18 at 10:31
1

The restriction on b.value = 'very_high' must be part of the join condition. As it is now, you are doing the spatial join, then you remove the rows whose b.value is not very_high but this value is null for polygons that do not intersect. By putting the restriction in the join condition, you still get the data from the g table with no matching b rows.

At last, you want to use coalesce to replace null result with 0.

select g.gid, g.name as subdistrict,   
       coalesce(round(ST_area(ST_union(ST_INTERSECTION(g.geom,b.geom)))::numeric/100/g.area_ha, 2),0) as "share_value_very_high",          
       ST_union(ST_INTERSECTION(g.geom,b.geom)) as geom_union
FROM structure.subdistrict as g
left join soil.soilmap AS b 
on (ST_Intersects(b.geom, g.geom)  AND b.value = 'very_high')
group by g.gid, g.name, coalesce(b.value,'other');
  • I see the mistake, now it works perfectly! Thank you! – oppia04 Sep 26 '18 at 13:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.