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

Try this approach: https://gis.stackexchange.com/questions/109692/how-to-replicate-arcgis-intersect-in-postgis — dr_jts, Dec 09 '20 at 17:03
And see this: http://blog.cleverelephant.ca/2019/07/postgis-overlays.html — dr_jts, Dec 09 '20 at 17:03

BERA · Accepted Answer · 2020-12-15T08:49:59.620

Since the number of classes are known and only eight 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 to pivot.

(But with some python pandas the pivoting gets alot easier:

import pandas as pd
import psycopg2
con = psycopg2.connect(database="lmv", user="postgres", password="somepassword", host="localhost")
sql = "select tk.ogc_fid plots, (st_area(st_intersection(ok.geom, tk.geom))/10000)::int area, kkod::int
       from ok_my_riks ok join tk_rutnat tk
       on st_intersects(ok.geom, tk.geom) where tk.ogc_fid in (5425, 5424, 5419, 5420)"
df = pd.read_sql_query(sql, con) #https://stackoverflow.com/questions/27884268/return-pandas-dataframe-from-postgresql-query-with-sqlalchemy
d = {601:'forest',303:'rural',901:'water',611:'grassland'} #From class code to class description
df['class'] = df.kkod.map(d) #https://stackoverflow.com/questions/29794959/pandas-add-new-column-to-dataframe-from-dictionary
df2 = df.groupby(by=['plots','class'])['area'].sum().reset_index().pivot(index='plots', columns='class', values='area').fillna(0).astype(int)

)

score 2 · Answer 2 · answered Dec 09 '20 at 17:08

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 Answers2