# Counting geometries which intersect another geometry using PostGIS

I'm trying to calculate the total area for each layer I have within a specific area and count the geometries which intersects with that specific area.

I can calculate the area and the results seem plausible, but how do I count the objects?

SELECT st_area(ST_Transform(ST_Setsrid((st_intersection(
st_union(st_makevalid(object.postgis_area_geometry)),
specific_area_geometry
)), 3857), 4326)::geography),
layer.name

FROM    object
LEFT JOIN layer ON layer.id = object.layer_id
WHERE
ST_ISVALID(object.postgis_area_geometry)
GROUP BY
layer.id;

for instance:

Here I have 4 Layers. Each layer has just one Geometry, just to keep it simple.

Object 1: Black 100 m2 = Layer 1
Object 2: Blue 70 m2 =   Layer 2
Object 3: Green 120 m2 = Layer 3
Object 4: Purple 70 m2 = Layer 4

The red square is the specific area within I want to know the total area of each layer.

So I want the following results:

Layer 1 = 35 m2, count = 1
Layer 2 = 10 m2, count = 1
Layer 3 = 35 m2, count = 1
Layer 4 = 0  m2, count = 0
• Any query with an ST_SetSRID is immediately suspect, especially in the context of an ST_Transform. You seem to have the IDs encoded already, so isn't this just a count(*) issue? Oct 22, 2021 at 12:11
• I just need to convert from EPSG 3857 to 4326 to get the area in square meter
– BR75
Oct 22, 2021 at 12:25
• If the data is already in 3857, then the SetSRID is a waste of computing resources. If you haven't set any SRID, then you should do that first (permanently, so it isn't necessary in subsequent queries) Oct 22, 2021 at 16:47
• Problem was that my postgis_area_geometry column is auto generated from another column I have and the default srid is 4326 which is wrong in my case because my data is in 3857, so this is because I had to use setSrid. I just fix this, thx.
– BR75
Oct 28, 2021 at 9:30
• Right. But we use a Question/Answer model, and answers don't belong in Questions. The edit can be rolled back, but others can't add an Answer for you. We don't do thanks here as a matter of policy see the Help Center. Oct 28, 2021 at 11:22

What you want is a combination of ST_Intersects and ST_Intersection. The count you're doing doesn't make a lot of sense to me, since you'd only ever get either 1 or 0. For that reason, you could just use ST_Intersects(a.geom,b.geom) for a boolean return of either True or False. And as per geozelot's comment, you can cast bool to int to get the output as you specify. I've also included how you'd get the full count of intersecting features.

I've created two layers; gisse (5 pink features) and gisse2 (1 blue feature), both using SRID: 4326:

select id,st_area(st_transform(geom,3857)) from gisse order by id

Gives:

id st_area
1 571990204493.9556
2 782180369787.4814
3 717504875654.8381
4 645453473213.9517
5 683918275314.2

To generate polygons showing the parts of table 1 that intersect with table 2:

create view theintersects as
select a.id,
st_area(st_transform(st_intersection(a.geom,b.geom),3857)),
st_intersection(a.geom,b.geom) from gisse a
inner join gisse2 b on st_intersects(a.geom,b.geom)
order by a.id

Computed features in yellow:

If I instead want to get the statistics without the geometry field, including non-intersecting rows:

select a.id,
st_area(st_transform(st_intersection(a.geom,b.geom),3857)),
st_intersects(a.geom,b.geom)::int as intersect_count
from gisse a left join gisse2 b on st_intersects(a.geom,b.geom)
order by a.id
id st_area intersect_count
1 0 0
2 56153016337.158646 1
3 261289632649.45654 1
4 49506223558.761505 1
5 197914797339.4397 1

If you just want the count of features from gisse that intersect with gisse2, execute:

select count(a.id) as intersecting_features
from gisse a inner join gisse2 b on st_intersects(a.geom,b.geom)

Result:

intersecting_features
4
• I just edit my question with a little more details of what I want to achieve.
– BR75
Oct 22, 2021 at 12:42
• Sure, I'll expand on the answer. What you want is a combination of ST_Intersects and ST_Intersection Oct 22, 2021 at 13:41
• FYI: PostgreSQL happily casts BOOL to INT and vice versa, i.e. SELECT ST_Intersects(a.geom, b.geom)::INT. Oct 26, 2021 at 14:37
• Thanks @geozelot - updated the answer with your suggestion. Oct 27, 2021 at 14:36
• Sorry for the confusion, I just want to count all objects, for each layer which intersects with that specific area and calculate the total area for each layer within the specific area. It does not make sense with only one object in a layer, but I have a lot of objects in each layer. I think I got this, see my edit. The purpose of this is that this specific area is a user input. So basically he wants to get a report about what's going on in an area he selects. In this case I show him the number of objects for each layer and the total area for each layer.
– BR75
Oct 28, 2021 at 10:18