PostGIS multiple spatial joins in the same query

Question

While I understand how I would do something like this in non-spatial SQL, I have confused about how I would join data from more than two joins in one query in spatial SQL.

I have a table of zipcode polygons 'Z', a table of restaurant polygons 'F' and restaurant points 'FP'. What I would like to produce as a result is a table with a row for each zip code and columns for the count of the restaurant polygons, the count of restaurant points, and the sum of both restaurant columns within each zip code.

As a test, I am only using 1 zipcode. First just getting the restaurant points:

SELECT z.modzcta,COUNT(f.name)
FROM osm.food_point f
INNER JOIN osm.modzcta z
ON ST_Intersects(ST_Transform(z.geom,3857),f.way)
WHERE z.modzcta = '11238'
GROUP BY z.modzcta;

returns [('11238', 192)]

Then, using the polygons:

SELECT z.modzcta,COUNT(fp.name)
FROM osm.food_polygon fp
INNER JOIN osm.modzcta z
ON ST_Intersects(ST_Transform(z.geom,3857),fp.way)
WHERE z.modzcta = '11238'
GROUP BY z.modzcta;

I get [('11238', 10)]

When I try to combine these two together:

SELECT z.modzcta,COUNT(fp.name) res_points, COUNT(f.name) res_polygons, COUNT(fp.name) + COUNT(f.name) total_rest
FROM osm.food_polygon f
INNER JOIN osm.modzcta z
ON ST_Intersects(ST_Transform(z.geom,3857),f.way)
INNER JOIN osm.food_point fp
ON ST_Intersects(ST_Transform(z.geom,3857),fp.way)
WHERE z.modzcta = '11238'
GROUP BY z.modzcta;

I get [('11238', 1920, 1930, 3850)]

I'm not sure what is going on here.

Also as an aside, is there any way to return the column names with the result of a query in the psycopg2 python module?

Answer 1

Your JOIN chain multiplies rows of the middle table for every pair of rows of the right-most table (i.e. 192 point * 10 polygons = 1920 total point counts).

While DISTINCT works in this case (it may not in others), the JOIN chain is still highly inefficient due to the unnecessary duplication; better use a correlated sub-query, or better even, a CROSS JOIN LATERAL:

SELECT z.modzcta, r_pnt.cnt, r_ply.cnt
FROM   osm.modzcta AS z
CROSS JOIN LATERAL (
  SELECT COUNT(f.*) AS cnt
  FROM   osm.food_point AS f
  WHERE  ST_Intersects(ST_Transform(z.geom,3857),f.way)
) AS r_pnt
CROSS JOIN LATERAL (
  SELECT COUNT(f.*) AS cnt
  FROM   osm.food_polygon AS f
  WHERE  ST_Intersects(ST_Transform(z.geom,3857),f.way)
) AS r_ply
;

Make sure there is a GIST index on the geometry columns of both join tables.