Count number of different points within all polygons?

Question

I have 2 PostGIS tables in PostgreSQL: polytable and pointtable. The point table has 1000s of points and a class column with 6 different unique classes (integer values from 1 to 6). The polygon table has 1000s of polys and a name column. Some names appear multiple times.

For every single polygon, I want to count how many points (per point class) are in each. Although the polygons contain a name column (we'll call polyname here), I don't need these grouped at all... The current query I use only works when using GROUP BY polyname on the polygon table. Therefore I can't treat each polygon individually:

SELECT
  polyname,
  count(pointclass) FILTER (WHERE pointclass = 1) AS pointclass1,
  count(pointclass) FILTER (WHERE pointclass = 2) AS pointclass2,
  count(pointclass) FILTER (WHERE pointclass = 3) AS pointclass3,
  count(pointclass) FILTER (WHERE pointclass = 4) AS pointclass4,
  count(pointclass) FILTER (WHERE pointclass = 5) AS pointclass5,
  count(pointclass) FILTER (WHERE pointclass = 6) AS pointclass6
FROM 
    (select poly.polyname, pnt.pointcolumn
    from schema.polytable poly
    left join schema.pointtable pnt
    on st_intersects(pnt.geom, poly.geom)) sub
GROUP BY polyname

How can I remove the grouping in the query above? Simply removing the last line throws an error (column "sub.polyname" must appear in the GROUP BY clause or be used in an aggregate function)

This is my desired output:

polyid   polyname   polygeom   pointclass1   pointclass2   pointclass3...
1        red        0101000..  55            10            5
2        blue       0101000..  3             25            46
3        green      0101000..  100           39            11
4        yellow     0101000..  0             0             3
5        red        0101000..  0             11            37
6        green      0101000..  58            99            200
7        green      0101000..  0             0             0

Note above, how I still want polygon rows returned which have no points within. So a LEFT JOIN of point count to the polygon, per point class. Similar to the query below, but with 6 point count columns for each point class:

SELECT poly.id, poly.name, count(point.shape) AS pointcount 
FROM polytable poly
LEFT JOIN pointtable point ON st_intersects(poly.geom, point.grom)
GROUP BY poly.id

I'd prefer not to take the longer route of adding 6 empty columns onto the polygon table for each point class, then running an UPDATE on these columns to populate them, using different WHERE pointclass = 1/2/3... etc... There must be a better way?

Answer 1

The solution is to group by both polyid and polyname, it should return all your desired records. Just call the polyid column in all your select statements including the subquery. You can also call the geometry of your polygons by adding the poly.geom.

SELECT
  polyid, polyname, poly.geom
  count(pointclass) FILTER (WHERE pointclass = 1) AS pointclass1,
  count(pointclass) FILTER (WHERE pointclass = 2) AS pointclass2,
  count(pointclass) FILTER (WHERE pointclass = 3) AS pointclass3,
  count(pointclass) FILTER (WHERE pointclass = 4) AS pointclass4,
  count(pointclass) FILTER (WHERE pointclass = 5) AS pointclass5,
  count(pointclass) FILTER (WHERE pointclass = 6) AS pointclass6
FROM 
    (select poly.polyid, poly.polyname, poly.polygeom, pnt.pointclass
    from schema.polytable poly
    left join schema.pointtable pnt
    on st_intersects(pnt.geom, poly.geom)) sub
GROUP BY polyid, polyname, poly.geom

Alternatively, you can use the following simpler query, since in this case, you don't even need to use a subquery:

select
  polyid, polyname, poly.geom,
   count(pointclass) FILTER (WHERE pointclass = 1) AS pointclass1,
   count(pointclass) FILTER (WHERE pointclass = 2) AS pointclass2,
   count(pointclass) FILTER (WHERE pointclass = 3) AS pointclass3,
   count(pointclass) FILTER (WHERE pointclass = 4) AS pointclass4,
   count(pointclass) FILTER (WHERE pointclass = 5) AS pointclass5,
   count(pointclass) FILTER (WHERE pointclass = 6) AS pointclass6
FROM schema.polytable poly
left join schema.pointtable pnt
on st_intersects(pnt.geom, poly.geom)
GROUP BY polyid, polyname, poly.geom