This is essentially a duplicate question of multiple others, with the sole difference being a table self-join.
However, all queries currently present in this post have delicate CRS misunderstandings, at least when it comes to distances:
- the main problem here is the threshold given to
ST_DWithin; the units of that value are CRS dependent, thus, as the data is in EPSG:4236, you are searching in a radius of 10000 degrees!
- follow-up problem is the actual distance a degree represents in suface distance; one degree of Longitude does not represent the same suface distance over different Latitudes!
Arguably the best way to realize true KNN searches uses the
LATERAL JOIN in conjunction with the
<-> KNN operator, and, optionally, a limiting radius filter (e.g.
&& BBox comparator).
Concerning the units, one could choose a on-the-fly cast to geography type to tackle the CRS/distance issues and get the most precise distances in one go, using speroidal (or, quicker, spherical) algebra.
ST_Distance(p1.geom::geography, p2.geom::geography) AS dist
FROM points AS p1
JOIN LATERAL (
WHERE ST_DWithin(p1.geom::geography, geom::geography, 10000)
AND p1.id <> id
ORDER BY p1.geom::geography <-> geom::geography
) AS p2
dist in meter to the nearest neighbor
p2.id for each
Note: due to the cast to geography, units for any accepting function will be in meter as well, thus the
As already mentioned, it is essential that you have a proper index in place! Checking the
EXPLAIN ANALYZE is crucial to find out if it is actually used (although you can tell it is if you get results within your lifetime I guess...), and running
VACUUM ANALYZE <table_name> in advance can help to enforce its use.
Now, the liberal use of the on-the-fly cast to geography will take a heavy toll on execution speed. I´d recommend to either project the data to a suitable projection for distance measurements of your area, or, possibly better, change the geometry type to geography; both can be achieved by adding a new column, if you want your original
geoms to stay untouched, and add an own index to it.
Using test data on 70.000 points with porperly indexed geography column (having added a second one) took about 1 min. to complete the initial, uncached run on a mid tech setup.