4

I am trying to extract information from a base map based on overlaying polygons as shown in the picture.

enter image description here

The base map represents different groups (Names 1-8).

From the overlaying polygons (gray,names 1-7) I would like to know the following:

What is the area or proportion of the overlay (e.g.Polygon3= 40% (purple) in group 1, 60% in group2 (orange))

I tried to get there using:

SELECT 
t1.name,
t2.group_name
FROM 
public.overlay AS t1,
public.basemap AS t2
WHERE st_intersects(t1.geom,t2.geom)=TRUE
ORDER BY t1.name;

This gives me the underlying groups.

However, as some geometries of the base map are in the same group (same color) it would be good to have this summed up (e.g. polygon 4 intersects three geometries of the base map which are all belonging to the purple group).

I did this using R before but would prefer to do it in the database directly.

My preferred output would look something like this :

name    |group_name| coverage
1       |1         | 1
2       |2         | 1
3       |1         | 0.4
3       |2         | 0.6
4       |1         | 1
5       |1         | 0.65
5       |2         | 0.35
...
2
  • 2
    It is GIS SE policy to only ask one question per Question. Having multiple questions, some of which are maybe answered in the question, is confusing, as is using custom functions without defining them (what does utmzone() do? Why would this be better than using a geography cast?). Please edit the question to focus on the question you want answered (Hint: the answer will involve ST_Intersection())
    – Vince
    Dec 29, 2016 at 19:17
  • @Vince Thank you for the comment. I edited the question and hope it is more clear now. Dec 29, 2016 at 20:50

1 Answer 1

3

Assuming that the SRID of the geometries is not geographic, this should do the trick:

SELECT t1.name,
       t2.group_name,
       sum(st_area(st_Intersection(t1.geom,t2,geom))/st_area(t1.geom)) as coverage
FROM    public.overlay t1
LEFT OUTER JOIN public.basemap t2 ON st_intersects(t1.geom,t2.geom)
GROUP BY t1.name,t2.group_name
ORDER BY t1.name;

If the geometries are in degrees, then you'll need to add a ::geography cast within the st_area functions to generate geodetic area (not Cartesian degrees).

Note that this would also handle any overlay which doesn't intersect any basemap elements.

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.