1

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?

5
  • 1
    That query tells you that a tiny fraction of your Lake is intersected by features from summ_table that have an area greater than the threshold. All others are, in sum, fully covered by features from summ_table that have an area greater than the threshold. Note that you are filtering rows from summ_table, not from the subset. Note also that you are also intersecting the same polygons with each other except for Lake.
    – geozelot
    Commented Nov 13, 2020 at 20:08
  • Are you trying to filter on the sums of intersected areas per subset county that are greater than 7769967522? If so, then maybe try a HAVING clause at the very end, instead of a WHERE clause, like this,: HAVING sum(inter.area) > 7769967522, otherwise your filter merely checks the total area of the counties, not the intersected areas.
    – FSimardGIS
    Commented Nov 14, 2020 at 3:03
  • @geozelot - thank you for this explanation. I am still confused as to how Lake (as a part of zones) could intersect with anything in summ_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 the ST_Intersects portion of the query with ON (ST_Relate(summ_table.geom, zones.geom, 'T********')) does not change the result
    – rumski20
    Commented Nov 17, 2020 at 15:04
  • @FSimardGIS - thank you for this suggestion. I will attempt the approach you suggested.
    – rumski20
    Commented Nov 17, 2020 at 15:05
  • Quite like you found out; it suggests there are tiny overlapping areas, probably due to topological mismatches, or floating point rounding errors. Also, as I said, with that query you are also getting the intersection between the same polygons (same id in both tables), so the covered area will always be 100%. Use summ_table.ogc_fid <> zone.ogc_fid in your WHERE filter to not consider the same geometries.
    – geozelot
    Commented Nov 17, 2020 at 16:37

1 Answer 1

0

After guidance from @geozelot and @FSimardGIS, I took another look at the geometry involved with this query.

select 
    row_number() OVER() As id, 
    st_setsrid(st_collectionextract(st_intersection(lake.geom, slc.geom), 1), 26915) as geom
from (
    select geom from mn_county_dissolved where id = 69
) as slc,
(
    select geom from mn_county_dissolved where id = 38
) as lake

enter image description here

Revealed that there is somehow some sliver geometry that is overlapping between these two featueres which explains the odd result from the query above.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.