Postgres/PostGIS spatial join returning null results

Question

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.

Answer 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)