7

In PostgreSQL/PostGIS I have 2 tables:

  • A POLYGON table (polytable1) with shapes named 'Poly 1', 'Poly 2' etc... Each name can appear multiple times (i.e. 20 polygons with the name 'Poly 1')

  • A POINT table (pointtable) where each point has a simple letter name ('w','x','y' etc) There are only 5 variations of name (only 4 shown in this example).

I want to create a copy of the polytable table (or just run as a SELECT output) but with new columns storing a count of how often a point falls inside it. The columns can be named after the point names.

To illustrate for an example, I'd like to take this scenario for polytable1 and pointtable:

postgis

and produce this output for polytable2/SELECT view:

enter image description here

  • Do you know all possible point letters? How many are they? – BERA Apr 2 at 16:56
  • @BERA yes there are actually just 5 letters in total. V,W,X,Y,Z – Theo F Apr 2 at 16:58
  • Since you get those errors; what is the SRID of both geom columns? This is not an issue of any of the queries here. – geozelot Apr 2 at 19:18
  • Both are EPSG 27700 – Theo F Apr 2 at 19:37
9

Try:

SELECT
  polyname,
  count(pid) FILTER (WHERE pid='w') AS "w",
  count(pid) FILTER (WHERE pid='x') AS "x",
  count(pid) FILTER (WHERE pid='y') AS "y",
  count(pid) FILTER (WHERE pid='z') AS "z"
FROM 
    (select polygons.polyname,points.pid
    from polygons
    left join points
    on st_intersects(points.geom, polygons.geom)) sub  
GROUP BY polyname

enter image description here enter image description here

| improve this answer | |
  • 1
    A slight variation on this using JOIN LATERAL is: SELECT polyname, count(pid) FILTER (WHERE pid='w') AS "w", count(pid) FILTER (WHERE pid='x') AS "x", count(pid) FILTER (WHERE pid='y') AS "y", count(pid) FILTER (WHERE pid='z') AS "z" FROM polygons LEFT JOIN LATERAL ( SELECT pid FROM points WHERE ST_Intersects(points.geom, polygons.geom) ) AS sub ON true GROUP BY polyname – dr_jts Apr 2 at 18:01
  • @dr_jts Would that be faster or why use lateral join? – BERA Apr 2 at 19:02
  • @BERA I get a ERROR: Operation on mixed SRID geometries CONTEXT: parallel worker SQL state: XX000 using your solution. – Theo F Apr 2 at 19:06
  • @BERA, hard to say about which is faster. Probably depends on the data. – dr_jts Apr 2 at 20:44
  • As for your query, could the inner SELECT be removed, like this: SELECT polyname, count(pid) FILTER (WHERE pid='w') AS "w", count(pid) FILTER (WHERE pid='x') AS "x", count(pid) FILTER (WHERE pid='y') AS "y", count(pid) FILTER (WHERE pid='z') AS "z" FROM polygons LEFT JOIN points ON st_intersects(points.geom, polygons.geom) GROUP BY polyname – dr_jts Apr 2 at 20:45
4

Not exactly what you want to get, the following query will create a table with polyname, pointname, count

SELECT polyname, pointname, count(*) as cnt 
FROM poly_table inner join point_table on (ST_Contains(poly_table.geom, point_table.geom)
GROUP BY polyname, pointname
ORDER BY polyname, pointname

A supposed the name of the geometry column in both table is "geom". This way the zero count won't be stored.

The result should be:

Poly1   w   3
Poly1   x   6
Poly1   z   1
Poly2   w   6
Poly2   y   2
Poly3   w   2
...
| improve this answer | |
4

I think dynamic column names would be quite difficult. But if they don't need to be dynamic it sounds pretty simple:

select
po.poly_name,
sum(case when pt.point_name = 'w' then 1 else 0 end) as w,
sum(case when pt.point_name = 'x' then 1 else 0 end) as x,
sum(case when pt.point_name = 'y' then 1 else 0 end) as y,
sum(case when pt.point_name = 'z' then 1 else 0 end) as z
from
polytable1 po
left join pointtable pt on st_intersects(po.geom, pt.geom)
group by po.poly_name

Edit for OP:

select
po.poly_name,
sum(case when pt.point_name = 'w' then 1 else 0 end) as w,
sum(case when pt.point_name = 'x' then 1 else 0 end) as x,
sum(case when pt.point_name = 'y' then 1 else 0 end) as y,
sum(case when pt.point_name = 'z' then 1 else 0 end) as z
from
polytable1 po
left join pointtable pt on st_intersects(st_setsrid(po.geom, 27700), st_setsrid(pt.geom, 27700))
group by po.poly_name
| improve this answer | |
  • column names will be static. – Theo F Apr 2 at 17:13
  • I also get a ERROR: Operation on mixed SRID geometries CONTEXT: parallel worker SQL state: XX000 error with this. – Theo F Apr 2 at 19:10
  • One of your tables is not 27700, but if you say it is, then try the updated query. – HeikkiVesanto Apr 2 at 20:47
  • ah it seems that one of my tables has been uploaded as SRS:900914. No idea where that came from! I'll retry – Theo F Apr 2 at 21:23
4

You can use the following code (I work in the schema public):

-- Creation of a new table : I use distinct to get one row for each value of polyname
CREATE TABLE polytable2 AS
SELECT DISTINCT polyname 
FROM polytable1 
ORDER BY polyname; 

-- Creation of four new fields : one for each letter with 0 as default value
ALTER TABLE polytable2 
ADD COLUMN w INT DEFAULT 0,
ADD COLUMN x INT DEFAULT 0,
ADD COLUMN y INT DEFAULT 0,
ADD COLUMN z INT DEFAULT 0;

-- Updating my 4 fields using a subquery that calculates the point count in each polyname value 
-- and then making a join between my two polygon tables. 
-- 4 updates, each time I change the name of the column and my letter in the where clause.
-- for the letter w
UPDATE polytable2
SET w = count FROM (SELECT polyname, count(*)
FROM polytable1, pointtable
WHERE letter = 'w' AND st_intersects(polytable1.geom, pointtable.geom)
GROUP BY polyname) as temp WHERE polytable2.polyname = temp.polyname ;

-- for the letter x
UPDATE polytable2
SET x = count FROM (SELECT polyname, count(*)
FROM polytable1, pointtable
WHERE letter = 'x' AND st_intersects(polytable1.geom, pointtable.geom)
GROUP BY polyname) as temp WHERE polytable2.polyname = temp.polyname;

-- for the letter y
UPDATE polytable2
SET y = count FROM (SELECT polyname, count(*)
FROM polytable1, pointtable
WHERE letter = 'y' AND st_intersects(polytable1.geom, pointtable.geom)
GROUP BY polyname) as temp WHERE polytable2.polyname = temp.polyname ;

-- for the letter z
UPDATE polytable2
SET z = count FROM (SELECT polyname, count(*)
FROM polytable1, pointtable
WHERE letter = 'z' AND st_intersects(polytable1.geom, pointtable.geom)
GROUP BY polyname) as temp WHERE polytable2.polyname = temp.polyname ;

I can explain more if you need.

| improve this answer | |

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