Postgres 13

By using db<>fiddle, you agree to license everything you submit by Creative Commons CC0.
-- original table
CREATE TABLE tbl (id int, lon numeric, lat numeric);
INSERT INTO tbl VALUES
  (1, 	10.111, 	20.415)  -- 1 neighbours at tol 0.0011
, (2, 	10.099, 	30.132)  -- 0 neighbours at tol 0.0011
, (3, 	10.110, 	20.414)  -- 1 neighbours at tol 0.0011
, (4, 	11.110, 	25.414)  -- 0 neighbours at tol 0.0011
, (11, 	12.111, 	20.415)  -- 1 neighbours at tol 0.0011
, (12, 	12.112, 	20.415)  -- 2 neighbours at tol 0.0011
, (13, 	12.113, 	20.414)  -- 1 neighbours at tol 0.0011
;

-- spatial espression index
CREATE INDEX tbl_point_gist_idx ON tbl USING GiST (point(lon, lat));
7 rows affected

-- basic nearest neighbour query
SELECT t.*, t1.id AS near_id, t1.lon AS near_lon, t1.lat AS near_lat
FROM   tbl t
JOIN   LATERAL (
   SELECT *
   FROM   tbl t1
   WHERE  t1.id <> t.id
   ORDER  BY point(t.lon, t.lat) <-> point(t1.lon, t1.lat)
   LIMIT  5  -- arbitrary
   ) t1 ON abs(t.lon - t1.lon) < 0.0011
       AND abs(t.lat - t1.lat) < 0.0011;


id
lon
lat
near_id
near_lon
near_lat


1
10.111
20.415
3
10.110
20.414

3
10.110
20.414
1
10.111
20.415

11
12.111
20.415
12
12.112
20.415

12
12.112
20.415
11
12.111
20.415

12
12.112
20.415
13
12.113
20.414

13
12.113
20.414
12
12.112
20.415


-- original task
WITH RECURSIVE d(tol) AS (SELECT float8 '0.0011')  -- input tolerance here
, cte AS (
   SELECT t.id, t.lon, t.lat, d.tol
        , 1 AS step, ARRAY [t1.id] AS neighbours
   FROM   tbl t
   CROSS  JOIN d
   JOIN   LATERAL (
      SELECT *
      FROM   tbl t1
      WHERE  t1.id <> t.id
      ORDER  BY point(t.lon, t.lat) <-> point(t1.lon, t1.lat)
      LIMIT  1
      ) t1 ON abs(t.lon - t1.lon) < d.tol
          AND abs(t.lat - t1.lat) < d.tol

   UNION ALL
   SELECT t.id, t.lon, t.lat, t.tol
        , t.step + 1, neighbours || t1.id
   FROM   cte t
   JOIN   LATERAL (
      SELECT *
      FROM   tbl t1
      WHERE  t1.id <> t.id
      ORDER  BY point(t.lon, t.lat) <-> point(t1.lon, t1.lat)  -- !
      OFFSET t.step                                            -- !
      LIMIT  1                                                 -- !
      ) t1 ON abs(t.lon - t1.lon) < t.tol                      -- !
          AND abs(t.lat - t1.lat) < t.tol                      -- !
   )
SELECT DISTINCT ON (id) id, lon, lat, neighbours
FROM   cte
ORDER  BY id, step DESC;


id
lon
lat
neighbours


1
10.111
20.415
{3}

3
10.110
20.414
{1}

11
12.111
20.415
{12}

12
12.112
20.415
{11,13}

13
12.113
20.414
{12}


-- better: check actual geometric distance
WITH RECURSIVE d(tol) AS (SELECT float8 '0.001415')  --  input tolerance here
, cte AS (
   SELECT t.id, t.lon, t.lat, d.tol
        , 1 AS step, ARRAY [t1.id] AS neighbours
   FROM   tbl t
   CROSS  JOIN d
   JOIN   LATERAL (
      SELECT *, point(t.lon, t.lat) <-> point(t1.lon, t1.lat) AS dist
      FROM   tbl t1
      WHERE  t1.id <> t.id
      ORDER  BY dist
      LIMIT  1
      ) t1 ON dist < d.tol

   UNION ALL
   SELECT t.id, t.lon, t.lat, t.tol
        , t.step + 1, neighbours || t1.id
   FROM   cte t
   JOIN   LATERAL (
      SELECT *, point(t.lon, t.lat) <-> point(t1.lon, t1.lat) AS dist
      FROM   tbl t1
      WHERE  t1.id <> t.id
      ORDER  BY dist         -- !
      OFFSET t.step          -- !
      LIMIT  1               -- !
      ) t1 ON dist < t.tol   -- !
   )
SELECT DISTINCT ON (id) id, lon, lat, neighbours
FROM   cte
ORDER  BY id, step DESC;


id
lon
lat
neighbours


1
10.111
20.415
{3}

3
10.110
20.414
{1}

11
12.111
20.415
{12}

12
12.112
20.415
{11,13}

13
12.113
20.414
{12}
id	lon	lat	near_id	near_lon	near_lat
1	10.111	20.415	3	10.110	20.414
3	10.110	20.414	1	10.111	20.415
11	12.111	20.415	12	12.112	20.415
12	12.112	20.415	11	12.111	20.415
12	12.112	20.415	13	12.113	20.414
13	12.113	20.414	12	12.112	20.415