I am running an intersect from a subset of one table on itself (this is part of a larger parameterized query but the issue in question arose out of this particular situation) using a CTE (zones
in the query below).
zones
should result in four records selected by ID.
In the subsequent query summ_table
should be limited to just three of those four records due to the where clause (WHERE (ST_Area(summ_table.geom) > 7769967522)
).
However, when I run the query it produces the following result:
id area
4 7911449533.18243
36 8167237869.65092
38 5.16730218436123e-05
69 17451037319.6371
id=38
is obviously not correct. I assumed the result would only include the three records that satisfy the above where clause.
Here is the full query:
-- select from zone table using IDs
WITH zones AS
(
SELECT ogc_fid, geom
FROM mn_county_dissolved
WHERE ogc_fid IN (4,36,69,38)
)
-- individual area of features that intersect each zone
SELECT inter.ogc_fid, sum(inter.area) AS area
FROM
(
-- zone / summary table intersect
SELECT
zones.ogc_fid,
CASE
WHEN ST_CoveredBy(summ_table.geom, zones.geom) THEN ST_Area(summ_table.geom)
ELSE ST_Area( ST_Intersection(summ_table.geom, zones.geom) )
END AS area
FROM mn_county_dissolved AS summ_table
INNER JOIN zones
ON ( ST_Intersects(summ_table.geom, zones.geom) )
WHERE (ST_Area(summ_table.geom) > 7769967522)
) AS inter
GROUP BY inter.ogc_fid
And here are the values for those records in question for ST_Area
:
ogc_fid | countyname | ST_Area
---------+-------------+------------------
4 | Beltrami | 7911449533.18243
36 | Koochiching | 8167237869.65092
38 | Lake | 5924950382.16346
69 | St. Louis | 17451037319.6371
I'll help you out with the math --> 7769967522 - 5924950382.16346 = 1845017139.84
As you can see the record with id=38
should definitely not satisfy that where clause
The only thing I can think of is that it might have something to do with the number of geometries for each record.
ST_NumGeometries
returns 8 for the record id=38
and 1 for all of the others.
Note that ST_GeometryType
returns ST_MultiPolygon
for all rows.
Why is that where clause in the query not filtering out that particular row (id=38
), and if the query is indeed working as expected, is it due to the fact that there is more than one geometry for that particular record?
Lake
is intersected by features fromsumm_table
that have an area greater than the threshold. All others are, in sum, fully covered by features fromsumm_table
that have an area greater than the threshold. Note that you are filtering rows fromsumm_table
, not from the subset. Note also that you are also intersecting the same polygons with each other except for Lake.HAVING sum(inter.area) > 7769967522
, otherwise your filter merely checks the total area of the counties, not the intersected areas.Lake
(as a part ofzones
) could intersect with anything insumm_table
when it is not a part of it. Should I be looking for some smaller sliver geometry? I think I can confirm that it is not a shared boundary issue since theST_Intersects
portion of the query withON (ST_Relate(summ_table.geom, zones.geom, 'T********'))
does not change the resultsumm_table.ogc_fid <> zone.ogc_fid
in yourWHERE
filter to not consider the same geometries.