I got 2 sets of points in 2 separate tables. Table_a got 100k points and table_b got 300k of points. I trying to find nearest points in relation find me any point from table_b that is within 50 meters from tabla_a. After that calculate fall column, group them by table_a a_id column and return highest value.

I wrote a following query that meet this criteira

         table_b.height - st_3ddistance(table_b.geom, table_a.geom) fall,
       FROM table_a
         INNER JOIN table_b ON _st_3ddwithin(table_a.geom, table_b.geom, 50)) a
WHERE fall >= 0
ORDER BY a_id, fall DESC;

I added 3d geometry indexes:

CREATE INDEX table_a_geom ON table_a USING GIST (geom gist_geometry_ops_nd);
CREATE INDEX table_b_geom ON table_b USING GIST (geom gist_geometry_ops_nd);

However my problem is that i can't make query to use them. Query planer is keep choosing sequence scan that is slow. I run some test changing _st_3ddwithin with st_3ddwithin, <<->> < 50 , creating 50 m buffer and intersect, st_3ddistance < 50 but everytime planner is choosing sequence scan. Is there a way to use indexes with higher performance or changing the query to use indexes?

My query plan:

Unique  (cost=10462593.70..10473018.43 rows=1 width=144)
  ->  Sort  (cost=10462593.70..10467806.06 rows=2084945 width=144)
        Sort Key: table_a.nmbayuid, ((table_b.height - st_3ddistance(table_b.geomgr, table_a.geom))) DESC
        ->  Nested Loop  (cost=0.00..10243762.28 rows=2084945 width=144)
              Join Filter: (_st_dwithin(table_a.geom, table_b.geomgr, '50'::double precision) AND ((table_b.height - st_3ddistance(table_b.geomgr, table_a.geom)) >= '0'::double precision))
              ->  Seq Scan on table_b  (cost=0.00..1459.47 rows=47147 width=96)
              ->  Materialize  (cost=0.00..10.97 rows=398 width=56)
                    ->  Seq Scan on table_a  (cost=0.00..8.98 rows=398 width=56)
  • 1
    What exactly is e_wires_mv12404 which is in the query plan but not the SQL? What does the query plan for just the inner query look like? I suggest not using function that start with _ST. Finally, you might be able to get better performance using ST_DWithin in 2D, using 35 meters, which is more or less the same as 50 meters from opposite edges of a cube. As you are looking for the single closest point within 50 meters, this might be a good candidate for a lateral join and using the ORDER BY a.geom <-> b.geom construct. Commented Jul 6, 2018 at 16:25
  • 1
    I had a similar problem last year, I dug up this post for you, let me know if it doesn't answer your questions?
    – WxGeo
    Commented Jul 6, 2018 at 17:38
  • 2
    If you look at the sql definition of the functions you see that the st_ functions like st_dwithin actually is a bounding box check and a call to the st function. It is the bounding box part that can use the index when you call the st function directly ther is no way for the database to use the index. You call the recheck function directly. Commented Jul 8, 2018 at 9:47
  • 1
    Would you like me to write up the lateral join solution, I think it would work well for what you describe Commented Jul 10, 2018 at 10:10
  • 1
    @AndreSilva functions starting with _ST are internal functions called by PostGIS after filtering with an index. If you call them directly, the index will not be used.
    – dbaston
    Commented Jul 10, 2018 at 17:21

3 Answers 3


Firstly, as has been noted in the comments, the leading underscore before ST function, ie, _ST_3DWithin will lead to the index not being used. I can't find any recent mention of this, but in older docs if you search for, eg, _ST_Intersects it states:

To avoid index use, use the function _ST_Intersects.

EDIT: As clarified by @dbaston in the comments, the functions with the leading underscore are internal functions that do not use the index when called and this continues to be the case (although it is hard to find in the docs).

Your query could possibly benefit from the LATERAL JOIN syntax, which lends itself well to k nearest neighbour (kNN) problems like this one.

   b.height - ST_3Ddistance(b.geom, a.geom) AS fall,
  FROM table_a a
          FROM table_b
          WHERE ST_3Ddwithin(a.geom, geom, 50)
          AND height - ST_3Ddistance(geom, a.geom) > 0
          ORDER BY height - ST_3Ddistance(b.geom, a.geom) DESC 
          LIMIT 1
        ) b ON TRUE;

This allows you to find the nearest k geometries from table a (in this case 1, due to LIMIT 1) to table b, ordered by the 3D distance between them. It is written using a LEFT JOIN, as it is conceivable that there might be some geometries in table a that are not within 50 meters of table b.

The lateral queries allow you to reference columns from the previous FROM clause, which makes it more powerful than standard sub queries, see the docs.

I can't test this against your data, but when I have run similar queries, the EXPLAIN statement indicates proper index use.

  • Your comments make great sens but I can't accept the answer because query you provided is doing different think that original query. As i mansion before " i am not looking for single closest point but a group of points within 50 meters and then i am selecting one with highest subtract value (height - ST_3Ddistance(geom, a.geom)) grouped by a_id
    – Losbaltica
    Commented Jul 11, 2018 at 7:00
  • I modified your query please have a look and add improvements if needed :)
    – Losbaltica
    Commented Jul 11, 2018 at 7:06
  • 1
    I modified the query, the only thing missing was "height -" in the order by. This will now find all points within 50 and return the one with the highest height - ST_3Ddistance(b.geom, a.geom) value. There is no need for distinct on, as this is all handled by each lateral query and LIMIT 1, ie, you will only get the largest fall value for each a_id. Commented Jul 11, 2018 at 8:00
  • Is this now working as you originally expected. Does EXPLAIN look sensible? Commented Jul 11, 2018 at 9:01
  • Is working as expected. Query performance is almost the same but the cost of the query is a lot smaller. New EXPLAIN: explain.depesz.com/s/Js5G I think i reach the limit of query optimisation and only think i can do now is to tune server or refactor the tables/logic. So it's answering me original question
    – Losbaltica
    Commented Jul 11, 2018 at 10:05

This link to PostGIS documentation recommends the following steps in order to ensure indexes and query planner are optimized:

  1. Make sure statistics are gathered about the number and distributions of values in a table, to provide the query planner with better information to make decisions around index usage. VACUUM ANALYZE will compute both.

  2. If vacuuming does not help, you can temporarily force the planner to use the index information by using the set enable_seqscan to off; command. This way you can check whether planner is at all capable to generate an index accelerated query plan for your query. You should only use this command only for debug: generally speaking, the planner knows better than you do about when to use indexes. Once you have run your query, do not forget to set ENABLE_SEQSCAN back on, so that other queries will utilize the planner as normal.

  3. If set enable_seqscan to off; helps your query to run, your Postgres is likely not tuned for your hardware. If you find the planner wrong about the cost of sequential vs index scans try reducing the value of random_page_cost in postgresql.conf or using set random_page_cost to 1.1;. Default value for the parameter is 4, try setting it to 1 (on SSD) or 2 (on fast magnetic disks). Decreasing the value makes the planner more inclined of using Index scans.

  4. If set enable_seqscan to off; does not help your query, it may happen you use a construction Postgres is not yet able to untangle. A subquery with inline select is one example - you need to rewrite it to the form planner can optimize, say, a LATERAL JOIN.

So, first try steps 1-3 before rewriting your query to use the indices. If that doesn't work, you could try to modify the query.

I believe (to the best of my ability to whip up SQL without running the code) that the query below will return identical results to yours, but don't know if it will be more efficient.

    table_b.b_id as b_id,
    table_b.height - st_3ddistance(table_b.geom, table_a.geom) as fall,
    table_b.geom as b_geom,
    table_a.a_id as a_id
    FROM table_a
         INNER JOIN table_b ON _st_3ddwithin(table_a.geom, table_b.geom, 50)) a
WHERE fall >= 0
ORDER BY a_id, fall DESC;
  • Very interesting after changing _st_3ddwithin to st_dwithin like other comments suggested and running VACUUM ANALYZE after, planner finally starts to catch the index !
    – Losbaltica
    Commented Jul 9, 2018 at 12:05

If you are using Postgres 10 (or newer), I would strongly recommend loading your data in Parallel tables.

You probably will need to spend time tuning it (data partitioning, and the number of workers), but I think is worth the effort. Theoretically, KNN is highly parallelizable, reaching constant time complexities, even O(1) if the amount of workers is equal to the number of elements where a KNN operation will be calculated.

Some practical reference on loading the data and performing the queries are provided here. He provides some detail on plan tunning (to force more workers to be actioned) here. Important to remark that parallel scripts involve a lot of task coordination, so that extreme theoretical bound of providing the most extreme parallelization doesn't hold in practice, due to networking, and other systems design characteristics.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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