# Counting features into intersections of Shapely Polygons

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: 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.

• Try this example. You can split the polygons for each 1-to-1 intersection and count each instance, make a list in python to populate with the count number and then the attribute table. – blu_sr Oct 10 '15 at 17:08
• Whilst there is no code included, but this answer on SO describes an algorithm to check if one polygon is entirely within another. I presume that if you checked for at least one line-intersection between line-segments then that would indicate overlapping polygons. – stevej Oct 12 '15 at 3:43
• Also there seems to be a geopandas function GeoSeries.intersects; I wonder if it works on polygons. – stevej Oct 12 '15 at 3:45
• do you have the ability to rasterize them? if you rasterize them all to have ones in the polygons, you can then use numpy to add them together and each number in the result will indicate how many polygons overlap at that pixel. – user1269942 Nov 29 '15 at 2:51

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);

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: 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
``````