2

I have a query where I'm choosing one polygon (parcel) that intersects a point. The problem I'm having is that when I run into a parcel that has multi-geometry, it often chooses the wrong polygon because I'm limiting the query to only return one result. How can I make this query "explode" this parcel and only choose the one polygon that intersects the point?

Here is my query:

SELECT  b.gid,
        st_astext(st_multi((st_dump(st_union(b.geom))).geom)) as area
FROM pointschema.point a
JOIN parcelschema.polygon b on st_within(a.geom, b.geom)
LIMIT 1

here is what the polygon looks like (you see the point in the NE corner): enter image description here

and here is what the result of the query is (orange): enter image description here

for reference I'm using:

  • PostgreSQL 9.5.1
  • PostGIS 2.2.1
  • Polygon geom: geom geometry(MultiPolygon,4326)
  • Point geom: geom geometry(Point,4326)
2
  • 1
    LIMIT 1 is unwise if you care about which feature is returned and don't provide an ORDER BY which will chose it. In the case of ST_Dump-ed geometries, you need to make the above query a virtual table (wrapped in a FROM ( ) alias) so you can reapply the search constraint on the parts.
    – Vince
    Nov 19, 2019 at 17:52
  • that's actually the path I am currently on - creating the result set in a CTE and basically duplicating it to only select from the dump.
    – Matt
    Nov 19, 2019 at 18:05

2 Answers 2

2

Here's a similar example of a point intersecting a 5-part multipolygon feature (from City of Denver Parcels - see below for SQL on finding one using ST_NumGeometries):

enter image description here

(maps from DBeaver spatial / geometry viewer)

First I create a CTE to dump out all the parts of the multipolygons:

select 
(ST_Dump(p.geom)).geom
    from ccd."CCD_Parcels" as p 
    where 1=1
and p.schednum = '0615201001999'

Then bring everything together:

;with cte_pt as (
    select st_setsrid(St_makepoint(-104.87635, 39.70721 ), 4326) as geom 
    )

, cte_poly as (

select 
(ST_Dump(p.geom)).geom
    from ccd."CCD_Parcels" as p 
    where 1=1
and p.schednum = '0615201001999'

)

select 
*
from
cte_pt as pt 
    join cte_poly as p 
on ST_within(pt.geom, p.geom)

Which gives me this result:

enter image description here

Now I'm not sure what exact result you're trying to get - are you just looking to return the part of the multipolygon feature that intersects the point and remove the rest?

I haven't used a LIMIT 1 or a ST_Union, but hoping this gives a bit more info on working with ST_Dump and multipolygons with intersects.

To find an example to work with, I used this SQL:

select * from (
select 
    schednum
, st_numgeometries(geom) as geomcount
from ccd."CCD_Parcels"
    ) as d 
   where d.geomcount = '5'
4
  • Thanks for your detailed answer, this will work great for my problem.
    – Matt
    Nov 19, 2019 at 19:17
  • @Matt also - I'm using DBeaver and the new spatial viewer for the maps... Nov 19, 2019 at 19:20
  • DBeaver looks amazing. I've been looking for a viewer where I can upgrade Postgres and not have to deal with the new PGadmin!
    – Matt
    Nov 19, 2019 at 19:24
  • Definitely try it out... ! It's come a long way! Nov 19, 2019 at 19:30
2

Since dumped geometries are not covered by the parent spatial index, here a slight improvement for large tables:

SELECT  py.gid,
        dmp.geom
FROM    (
    SELECT  a.gid,
            a.geom
    FROM    polygon AS a
    JOIN    point AS b
      ON    a.geom && b.geom
) AS    py,
        LATERAL ST_Dump(py.geom) AS dmp
JOIN    point AS pt
  ON    ST_Within(pt.geom, dmp.geom)
;

Here, the sub-select filters the polygon table using the index, to dump and relate only those polygons whose BBOX intersects the point(s BBOX)

2
  • Is this the same as doing an ST_Intersects with the points on the multipolygon data, then with that subset doing the same ST_Dump and ST_Within? Nov 19, 2019 at 21:54
  • @DPSSpatial yes that's the idea exactly. I used && on purpose, though, instead of ST_Intersects, as the latter would just add overhead.
    – geozelot
    Nov 20, 2019 at 7:34

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.