Filter is applied after DISTINCT ON

Question

There's a script at the end that will fully create the schema and populate it with sample data.

Schema

Consider these two tables:

Table of polygons:

CREATE TABLE my_polygon (
    my_polygon_id SERIAL PRIMARY KEY,
    common_id INTEGER NOT NULL,
    value1 NUMERIC NOT NULL,
    value2 NUMERIC NOT NULL,
    value3 NUMERIC NOT NULL,
    geom GEOMETRY(Polygon) NOT NULL
)
;

CREATE INDEX ON my_polygon (common_id);
CREATE INDEX ON my_polygon USING GIST (common_id, geom);

Table of points contained inside polygons:

CREATE TABLE my_point (
    my_point_id SERIAL PRIMARY KEY,
    common_id INTEGER NOT NULL,
    pointvalue NUMERIC NOT NULL,
    geom GEOMETRY(Point) NOT NULL
);

CREATE INDEX ON my_point (common_id);
CREATE INDEX ON my_point USING GIST (common_id, geom);

The fact that I'm using geometries isn't strictly related to the issue here; however, I think this makes the motives for what I'm attempting to do much clearer.

Problem query

The problem is that there are very tiny, insignificant overlaps between the polygons. (Trying to clean them is really not an option. The overlaps come from some kind of floating point error in generating them, as near as I can figure.) But some points can fall inside these tiny overlaps, resulting in two rows when I JOIN them based on containment. But really, each point should only be associated with a single polygon. When one does fall within two of them, it doesn't really matter which one it ends up associated with, so it's fine to make the query just pick one, like this:

SELECT DISTINCT ON (my_point.my_point_id)
    my_polygon.*,
    my_point.my_point_id,
    my_point.pointvalue,
    my_point.geom AS pointgeom
FROM my_polygon
JOIN my_point ON my_point.common_id = my_polygon.common_id AND ST_Contains(my_polygon.geom, my_point.geom)
WHERE my_polygon.common_id = 1
ORDER BY my_point.my_point_id, my_polygon.my_polygon_id

As in the above query, I normally want to SELECT based on the common_id. This query performs fine. It's query plan looks like this:

However, this is logic I need in a number of different queries, so I wanted to put it in a view. The result is that, as far as the query planner is concerned, the query looks like this:

SELECT *
FROM (
    SELECT DISTINCT ON (my_point.my_point_id)
        my_polygon.*,
        my_point.my_point_id,
        my_point.pointvalue,
        my_point.geom AS pointgeom
    FROM my_polygon
    JOIN my_point ON my_point.common_id = my_polygon.common_id AND ST_Contains(my_polygon.geom, my_point.geom)
    ORDER BY my_point.my_point_id, my_polygon.my_polygon_id
) point_with_polygon
WHERE common_id = 1

The result is that now PostgreSQL filters by common_id after performing the DISTINCT ON, which means that it has to JOIN the entirety of both tables. Here's it's query plan:

How can I allow PostgreSQL to push the filter down into an early part of the query and still put the general query in a view?

I'm stuck on PG 9.3 right now, but upgrading to 9.5 may be an option.

Schema and sample data script

Requires PostGIS. (That's why there's no SQL Fiddle.)

CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS btree_gist;

-- DROP FUNCTION ST_GeneratePoints(geometry, numeric);
DO $doblock$
BEGIN
    IF NOT EXISTS(SELECT * FROM pg_proc WHERE UPPER(proname) = UPPER('ST_GeneratePoints')) THEN
        -- Create naive ST_GeneratePoints if version of PostGIS is not new enough
        CREATE FUNCTION ST_GeneratePoints(g geometry, npoints numeric)
            RETURNS geometry
            VOLATILE
            RETURNS NULL ON NULL INPUT
            LANGUAGE plpgsql
            AS $$
            DECLARE
                num_to_generate INTEGER := npoints::INTEGER;
                adjustment CONSTANT FLOAT := 0.00000000001;
                x_min FLOAT := ST_XMin(g) + adjustment;
                x_max FLOAT := ST_XMax(g) - adjustment;
                y_min FLOAT := ST_YMin(g) + adjustment;
                y_max FLOAT := ST_YMax(g) - adjustment;
                temp_result GEOMETRY[];
                result_array GEOMETRY[] := ARRAY[]::GEOMETRY[];
            BEGIN
                IF ST_IsEmpty(g) THEN
                    RAISE EXCEPTION 'Cannot generate points inside an empty geometry';
                END IF;

                IF ST_Dimension(g) < 2 THEN
                    RAISE EXCEPTION 'Only polygons supported';
                END IF;

                -- Reduce number of loops to reduce slow array_cat calls
                WHILE num_to_generate > 0 LOOP
                    SELECT ARRAY_AGG(contained.point) INTO temp_result
                    FROM (
                        SELECT point
                        FROM (
                            SELECT ST_MakePoint(
                                x_min + random() * (x_max - x_min),
                                y_min + random() * (y_max - y_min)
                            ) point
                            -- Generate extras to reduce number of loops
                            --
                            -- Each point has a probability of  ST_Area(g) / ST_Area(ST_Envelope(g))  to fall within the polygon.
                            -- So on average, we expect ST_Area(g) / ST_Area(ST_Envelope(g)) of the points generated to fall within.
                            -- Generating  ST_Area(ST_Envelope(g)) / ST_Area(g) * num_to_generate  points means that on average, we'll
                            -- get
                            --
                            --     ST_Area(g) / ST_Area(ST_Envelope(g)) * ST_Area(ST_Envelope(g)) / ST_Area(g) * num_to_generate
                            --      = num_to_generate
                            --
                            -- points within the polygon. (Notice the numerators and denominators cancel out.) This means we'll 
                            -- only run one loop about half the time without generating an excessive number of points.
                            --
                            -- Generate at least 20 to avoid a lot of loops for small numbers, though.
                            FROM generate_series(1, GREATEST(20, CEIL(ST_Area(ST_Envelope(g)) / ST_Area(g) * num_to_generate)::INTEGER))
                        ) candidate
                        WHERE ST_Contains(g, candidate.point)
                        -- Filter out extras if we have too many matches
                        LIMIT num_to_generate
                    ) contained
                    ;
                    IF ARRAY_LENGTH(temp_result, 1) > 0 THEN
                        result_array := array_cat(result_array, temp_result);
                        num_to_generate := npoints - COALESCE(ARRAY_LENGTH(result_array, 1), 0);
                    END IF;
                END LOOP;
                RETURN (SELECT ST_Union(point) FROM UNNEST(result_array) result (point));
            END;
            $$;
        RAISE NOTICE 'Created ST_GeneratePoints';
    ELSE
        RAISE NOTICE 'ST_GeneratePoints exists';
    END IF;
END
$doblock$
;

DROP TABLE IF EXISTS my_polygon;

CREATE TABLE my_polygon (
    my_polygon_id SERIAL PRIMARY KEY,
    common_id INTEGER NOT NULL,
    value1 NUMERIC NOT NULL,
    value2 NUMERIC NOT NULL,
    value3 NUMERIC NOT NULL,
    geom GEOMETRY(Polygon) NOT NULL
)
;

CREATE INDEX ON my_polygon (common_id);
CREATE INDEX ON my_polygon USING GIST (common_id, geom);


WITH common AS (
    SELECT
        common_id,
        random() * 5000 AS common_x_translate,
        random() * 5000 AS common_y_translate
    FROM (
        SELECT TRUNC(random() * 1000) + 1 AS common_id
        FROM generate_series(1, 100)
        UNION 
        SELECT 1
    ) a
),
geom_set_with_small_overlaps AS (
    SELECT
        ST_MakeEnvelope(
            x.translate, 
            y.translate, 
            x.translate + 1.1, 
            y.translate + 1.1
        ) AS geom
    FROM
        generate_series(0, 9) x (translate),
        generate_series(0, 9) y (translate)
)
INSERT INTO my_polygon (common_id, value1, value2, value3, geom)
SELECT
    common_id,
    random() * 100,
    random() * 100,
    random() * 100,
    ST_Translate(geom, common_x_translate, common_y_translate)
FROM common, geom_set_with_small_overlaps
;

DROP TABLE IF EXISTS my_point;

CREATE TABLE my_point (
    my_point_id SERIAL PRIMARY KEY,
    common_id INTEGER NOT NULL,
    pointvalue NUMERIC NOT NULL,
    geom GEOMETRY(Point) NOT NULL
);

INSERT INTO my_point (common_id, pointvalue, geom)
SELECT
    common_id,
    random() * 100,
    (ST_Dump(ST_GeneratePoints(extent, FLOOR(5000 + random() * 15000)::NUMERIC))).geom
FROM (
    SELECT
        common_id,
        -- Small negative buffer prevents lying on the outer edge
        ST_Buffer(ST_Extent(geom), - 0.0001) AS extent
    FROM my_polygon
    GROUP BY common_id
) common
UNION ALL
SELECT
    common_id,
    random() * 100,
    (ST_Dump(ST_GeneratePoints(intersection, TRUNC(random() * 5)::NUMERIC))).geom
FROM (
    SELECT
        p1.common_id,
        p1.my_polygon_id AS id1,
        p2.my_polygon_id AS id2,
        ST_Intersection(p1.geom, p2.geom) AS intersection
    FROM my_polygon p1
    JOIN my_polygon p2 ON (
        p1.my_polygon_id < p2.my_polygon_id AND
        p1.common_id = p2.common_id AND
        ST_Intersects(p1.geom, p2.geom)
    )
) a
;

CREATE INDEX ON my_point (common_id);
CREATE INDEX ON my_point USING GIST (common_id, geom);

You probably want to VACUUM ANALYZE after that.

Query plans as text

WHERE clause inside (good performance):

Unique  (cost=1195.74..1207.74 rows=2400 width=216)
  ->  Sort  (cost=1195.74..1201.74 rows=2400 width=216)
        Sort Key: my_point.my_point_id, my_polygon.my_polygon_id
        ->  Nested Loop  (cost=5.34..1060.99 rows=2400 width=216)
              ->  Bitmap Heap Scan on my_polygon  (cost=4.93..191.74 rows=100 width=164)
                    Recheck Cond: (common_id = 1)
                    ->  Bitmap Index Scan on my_polygon_common_id_geom_idx  (cost=0.00..4.90 rows=100 width=0)
                          Index Cond: (common_id = 1)
              ->  Index Scan using my_point_common_id_geom_idx on my_point  (cost=0.41..8.68 rows=1 width=52)
                    Index Cond: ((common_id = 1) AND (my_polygon.geom && geom))
                    Filter: _st_contains(my_polygon.geom, geom)

WHERE clause outside (bad performance):

Subquery Scan on a  (cost=209447.85..215842.18 rows=1827 width=212)
  Filter: (a.common_id = 1)
  ->  Unique  (cost=209447.85..211274.80 rows=365390 width=212)
        ->  Sort  (cost=209447.85..210361.33 rows=365390 width=212)
              Sort Key: my_point.my_point_id, my_polygon.my_polygon_id
              ->  Nested Loop  (cost=0.41..63285.00 rows=365390 width=212)
                    ->  Seq Scan on my_polygon  (cost=0.00..338.00 rows=9800 width=164)
                    ->  Index Scan using my_point_common_id_geom_idx on my_point  (cost=0.41..6.41 rows=1 width=52)
                          Index Cond: ((common_id = my_polygon.common_id) AND (my_polygon.geom && geom))
                          Filter: _st_contains(my_polygon.geom, geom)

Answer 1

It turns out there's a way to work around this and get PG to optimize correctly: you have to include common_id in the DISTINCT ON clause.

Like this:

SELECT *
FROM (
    SELECT DISTINCT ON (my_point.my_point_id, my_polygon.common_id)
        my_polygon.*,
        my_point.my_point_id,
        my_point.pointvalue,
        my_point.geom AS pointgeom
    FROM my_polygon
    JOIN my_point ON my_point.common_id = my_polygon.common_id AND ST_Contains(my_polygon.geom, my_point.geom)
    ORDER BY my_point.my_point_id, my_polygon.common_id, my_polygon.my_polygon_id
) point_with_polygon
WHERE common_id = 1

This results in the same query plan as including the WHERE clause before the DISTINCT ON:

Unique  (cost=2307.77..2345.55 rows=5038 width=212)
  ->  Sort  (cost=2307.77..2320.36 rows=5038 width=212)
        Sort Key: my_point.my_point_id, my_polygon.my_polygon_id
        ->  Nested Loop  (cost=9.36..1479.97 rows=5038 width=212)
              ->  Bitmap Heap Scan on my_polygon  (cost=4.93..190.19 rows=100 width=164)
                    Recheck Cond: (common_id = 1)
                    ->  Bitmap Index Scan on my_polygon_common_id_geom_idx  (cost=0.00..4.90 rows=100 width=0)
                          Index Cond: (common_id = 1)
              ->  Bitmap Heap Scan on my_point  (cost=4.43..12.89 rows=1 width=52)
                    Recheck Cond: ((common_id = 1) AND (my_polygon.geom && geom))
                    Filter: _st_contains(my_polygon.geom, geom)
                    ->  Bitmap Index Scan on my_point_common_id_geom_idx  (cost=0.00..4.43 rows=2 width=0)
                          Index Cond: ((common_id = 1) AND (my_polygon.geom && geom))

Including common_id is kind of redundant since it's part of the JOIN condition anyway, but that fact also means that this won't change the result of the query.

Warning: Make sure everything matches up

It's very important that you use the common_id that's in the SELECT result. Take this query for example:

SELECT *
FROM (
    SELECT DISTINCT ON (my_point.my_point_id, my_point.common_id)
        my_polygon.*,
        my_point.my_point_id,
        my_point.pointvalue,
        my_point.geom AS pointgeom
    FROM my_polygon
    JOIN my_point ON my_point.common_id = my_polygon.common_id AND ST_Contains(my_polygon.geom, my_point.geom)
    ORDER BY my_point.my_point_id, my_point.common_id, my_polygon.my_polygon_id
) point_with_polygon
WHERE common_id = 1

This query uses my_polygon.common_id in the SELECT clause, but it uses my_point.common_id in the ORDER BY and DISTINCT ON clauses. PG will not push the filter down into the subquery in this case.