[GIS] Counting features into intersections of Shapely Polygons

geopandaspandaspythonshapely

I have a geopandas GeoDataFrame containing hundreds of shapely Polygon and MultiPolygon geometries. The polygons overlap in lots of places. I'd like to make a new geometry that contains a count of how many of them overlap. Something like this:

Counting overlaps

Does anyone have any ideas for how to approach this? I can't even see a way in.

Eventually I'd especially like to be able to weight the polygons, so that some polygons might be worth 2 on their own. Doing this with shapely's Z field might be nice.

Aside: I'm not particularly tied to any of these libraries, it's just where I've ended up. The coordinates in these geometries are actually pixel coordinates — I'm stumbling towards making a raster to overlay on another image. I'd prefer to keep my footprint as small as possible, since I'd like to be able to deploy this stuff to cloud servers etc, where I might not be able to install random stuff.

Best Answer

May be off topic because it is a postgresql/postgis solution:

In postgres/postgis it is a simple O(N^2) query that may/could be adopted to geopanda. 

$ psql gis_se;
-- create the postgis extension  
create extension postgis;

-- create a polygon table 
create table test_overlap(id serial primary key);

-- Add the geometry
select addgeometrycolumn('public','test_overlap','geom',4326,'POLYGON',2);
insert into test_overlap(geom) values
(ST_GeomFromEWKT('SRID=4326;POLYGON((-2 -2, -2 -1,-1 -1,-1 -2, -2 -2))')),
(ST_GeomFromEWKT('SRID=4326;POLYGON((-2 -2, -2 0, 0 0, 0 -2, -2 -2))')),
(ST_GeomFromEWKT('SRID=4326;POLYGON((2 2, 2 0, 0 0, 0 2, 2 2))')),
(ST_GeomFromEWKT('SRID=4326;POLYGON((2 2, 2 1,1 1,1 2, 2 2))')),
(ST_GeomFromEWKT('SRID=4326;POLYGON((-1.5 -1.5, -1.5 1.5,1.5 1.5,1.5 -1.5, -1.5 -1.5))'));

and defines 5 rectangles:

enter image description here

The intersection request with the table itself: 

select a.id, b.id, st_intersects(a.geom, b.geom) 
from test_overlap as a, test_overlap as b 
where a.id<>b.id;

shows which areas intersecting each other: 

 id | id | st_intersects 
----+----+---------------
  1 |  2 | t
  1 |  3 | f
  1 |  4 | f
  1 |  5 | t
  2 |  1 | t
  2 |  3 | t
  2 |  4 | f
  2 |  5 | t
  3 |  1 | f
  3 |  2 | t
  3 |  4 | t
  3 |  5 | t
  4 |  1 | f
  4 |  2 | f
  4 |  3 | t
  4 |  5 | t
  5 |  1 | t
  5 |  2 | t
  5 |  3 | t
  5 |  4 | t

Using this basis you can aggregate the counts for each ID object via the group by clausel:

select id, count(id)                         
from (select 
       a.id as id, b.id as bid, 
       st_intersects(a.geom, b.geom) as intersects 
       from test_overlap as a, test_overlap as b where a.id<>b.id
) as i
where intersects
group by id
order by id;

The result shows the wanted pattern. 

 id | count 
----+-------
  1 |     2
  2 |     3
  3 |     3
  4 |     2
  5 |     4
Related Question