How to get the area of two interesecting polygons divided into attributes on PostGIS

Question

I am trying to calculate the area overlap two vector polygon layers.
I have the one based layer - soil map and I have second layer - plots.
The example below:

I would like to know how much area each class takes.
I created a few queries (one for each class) and I merged to my statistics table. But this solution is too time consuming.
create table class_I as (
select e.id, sum(st_area(st_intersection(g.geom,e.geom))) as area 
from plots e, soils g 
where st_intersects(g.geom,e.geom) and g.class='I'
group by e.go_id);

I am looking for a solution that allows me to calculate all classes at once.
I would like my output to look something like this:
plots || area_classI || area_classII || area_classIII
  1   ||     0       ||   52864,28   ||     0
  2   ||    128      ||      0       ||  258687,89

BERA · Answer

Since you only have eight classes something like this is possible:
select tk.ogc_fid plots, 
            sum(st_area(st_intersection(ok.geom, tk.geom))) AS "total area",
            coalesce(sum(st_area(st_intersection(ok.geom, tk.geom))) FILTER (WHERE kkod=601), 0) AS "Forest", #kkod is the area class
            coalesce(sum(st_area(st_intersection(ok.geom, tk.geom))) FILTER (WHERE kkod=303), 0) AS "Rural",
            coalesce(sum(st_area(st_intersection(ok.geom, tk.geom))) FILTER (WHERE kkod=901), 0) AS "Water",
            coalesce(sum(st_area(st_intersection(ok.geom, tk.geom))) FILTER (WHERE kkod=611), 0) AS "Grassland"
            #more classes goes here

from ok_my_riks ok #My land classes table
    join tk_rutnat tk #My plot table
    on st_intersects(ok.geom, tk.geom)
    where tk.ogc_fid in (5425, 5424, 5419, 5420) #To limit the number of plots in my test
    group by plots

If unknown number of classes or high number of classes then I would try tablefunc (and quickly give up and read the select query into a pandas dataframe and pivot there instead)

dr_jts · Answer

This sounds like you want to compute the polygon coverage overlay of the two datasets.  The standard approach for doing this in PostGIS currently is to:

node the polygon linework using ST_Node or ST_Union (which is better to use in the latest release)
polygonize the noded lines using ST_Polygonize
determine resultant polygon parentage (if none indicates a hole) using ST_PointOnSurface and ST_Contains

This is well-described here.

How to get the area of two interesecting polygons divided into attributes on PostGIS

2 Answers

Add your own answers!

Ask a Question