2

I have 1.3 million points and 15 thousand polygons.

I'm trying to do a spatial join so I can find out how many points intersect with each polygon. Ultimately, I'll add up the values in some columns in the point data by each polygon.

So the point data is PCD4 and the polygon data is HCBG. I've verified that both tables have the same SRID. I tried this fairly simple spatial join query but the results are a null set which makes no sense!

select geoid, id
from public."PCD4"
join public."HCBG"
ON ST_Intersects(public."HCBG".geom, public."PCD4".geom)
group by id, public."PCD4".geom, geoid`

I've tried limiting the query to a single polygon but that also nulled out. I've checked the data in QGIS and there are definitely points in the polygons and specifically in the single test polygon I tried.

If I try to do the spatial join in QGIS, it errors out.

I've tested the integrity of the polygons. They're Census Tiger files and they checked out fine.

So I'm at a loss here. I'm sure the problem is something to do with the alignment of the data, either projection or something. But I am not sure how to troubleshoot it beyond checked for the SRID.

  • if you remove the public."PCD4".geom in the group by clause what happens? also if you had different SRID postgres would most likely tell you there are mixed SRIDs – ziggy Oct 25 '17 at 16:26
  • If I remove the geom from the group by I get the same result. – Matt Dempsey Oct 25 '17 at 16:28
  • try this...ST_Intersects(st_makevalid(public."HCBG".geom), st_makevalid(public."PCD4".geom)) – ziggy Oct 25 '17 at 16:30
  • 1
    ST_MakeValid shouldn't do anything in this case. I've run ST_IsValid on all the polygons and the points. There are no 'f' values for the returned query. So all my shapes and points are valid as is. – Matt Dempsey Oct 25 '17 at 16:38
  • 1
    seems to work for me... which version of postgres / postgis? Also, do the layers overlap visually when loaded into QGIS? This is to rule out getting the x and y values wrong way round when importing the point layer. If I had a penny for every time I'd done that... ;-) Probably not an issue if this is contiguous USA but worth checking. – Steven Kay Oct 25 '17 at 18:02
1

The more I use SQL, the more I find it simpler to break these queries down into their components, ie subqueries, instead of trying to do it all at once.
Let's start with a join to tag each point with the polygon it falls within. I'm using the st_covers function for this as it takes into account points that fall on a polygon boundary. I assign an alias to each table name to make them easier to reference in the query. I am not sure which ID goes with which layer, so I will make up names. Simply substitute with the correct ones.

SELECT pcd4.fid,
       hcbg.gid,
       pcd4.geom
FROM public."PCD4" pcd4
JOIN public."HCBG" hcbg
  ON st_covers(hcbg.geom, pcd4.geom)

This selects the ID of each point, and assigns the ID of the containing polygon. This can now be used as a WITH query for another SELECT query.

WITH pcd4_poly AS (
  SELECT pcd4.fid,
         hcbg.gid
  FROM public."PCD4" pcd4
  JOIN public."HCBG" hcbg
    ON st_covers(hcbg.geom, pcd4.geom)
)
SELECT hcbg.gid,
       count(pcd4_poly.fid),
       hcbg.geom
FROM public."HCBG" hcbg
WHERE hcbg.gid = pcd4_poly.gid
GROUP BY hcbg.gid

This query selects the ID and geometry of each polygon. The WITH command runs that select query to select the ID of each point, tagged with the ID of each polygon. When you select these together, it creates a one to many selection of polygons to points. Grouping on the polygon ID with a count of the point ID gives you the summary of point records per polygon that you are looking for.

I admit this may not be the most efficient way of doing it. It does, however, give you a process where it is easy to test each component to ensure they work separately, before putting them together.

To address what was at issue with your query, I notice you include what I think are the IDs from both layers, plus the point layer geometry in your GROUP BY clause. Since you are doing an intersection, there is no need to have a group clause as you aren't performing any aggregate operation in the SELECT query. The first part of your query would suffice. You may have had a NULL result simply because it was trying to display 1.3 million points.

select geoid, id
from public."PCD4"
join public."HCBG"
ON ST_Intersects(public."HCBG".geom, public."PCD4".geom)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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