計算地球上兩個經緯度位置之間的距離並不容易:要知道如何把地球當成球體處理,還要懂座標所屬的投影系統。所幸 contrib 內建的 earthdistance 擴充套件讓 PostgreSQL 輕鬆解決這個問題。

earthdistance contrib 擴充#

這類數學複雜到重新實作很容易出錯,所以我們要找已經被測試過的現成實作——earthdistance 正是為此而生:

create extension cube;
create extension earthdistance;

有了這個擴充,就能用它的 <@> 運算子,以 (longitude, latitude) 表示的 point 計算地表距離(單位是英里):

  select id, name, pos,
         round((pos <@> point(-0.12,51.516))::numeric, 3) as miles
    from pubnames
order by pos <-> point(-0.12,51.516)
   limit 10;
     id     │          name          │           pos           │ miles
════════════╪════════════════════════╪═════════════════════════╪═══════
   21593238 │ All Bar One            │ (-0.1192746,51.5163499) │ 0.039
   26848690 │ The Shakespeare's Head │ (-0.1194731,51.5167871) │ 0.059
  371049718 │ The Newton Arms        │ (-0.1209811,51.5163032) │ 0.047
  438488621 │ Marquis Cornwallis     │ (-0.1199612,51.5146691) │ 0.092
   21593236 │ Ship Tavern            │ (-0.1192378,51.5172525) │ 0.093
  ...
(10 rows)

最近的酒吧是 All Bar One,距離 0.039 英里(約 68.64 碼)。加上英里距離的計算幾乎沒有增加查詢時間——資料在記憶體中時,筆電上仍遠低於一毫秒。

酒吧與城市#

找最近的酒吧多容易,找最遠的就多容易——只要 order by ... desc。這個查詢本身未必有用,但它證明 kNN 搜尋也支援 ORDER BY DESC 的變體:

  select name, round((pos <@> point(-0.12,51.516))::numeric, 3) as miles
    from pubnames
order by pos <-> point(-0.12,51.516) desc
   limit 5;

接著想知道那些酒吧在哪個城市。用 OSM API 的以下 URL 可下載同一區域的城市清單:

http://www.overpass-api.de/api/xapi?*[place=city][bbox=-10.5,49.78,1.78,59]

微調 pubnames 專案的解析與匯入程式後,0.087 秒就匯入完城市名稱與位置:

create table if not exists cities
 (
   id   bigint,
   pos  point,
   name text
 );

create index on cities using gist(pos);

現在看看那些遙遠的酒吧位於何處——對每間酒吧找最近的已知城市,這種寫法叫關聯子查詢(correlated subquery)

   select name,
         (select name from cities c order by c.pos <-> p.pos limit 1) as city,
         round((pos <@> point(-0.12,51.516))::numeric, 3) as miles
    from pubnames p
order by pos <-> point(-0.12,51.516) desc
   limit 5;
      name       │  city  │  miles
═════════════════╪════════╪═════════
 Tig Bhric       │ Galway │ 440.194
 TP's            │ Galway │ 439.779
 Begley's        │ Galway │ 439.752
 Ventry Inn      │ Galway │ 438.962
 Fisherman's Bar │ Cork   │ 439.153
(5 rows)

同樣的查詢也能寫成 LATERAL join 結構,結果相同:

  select c.name as city, p.name,
         round((pos <@> point(-0.12,51.516))::numeric, 3) as miles
    from pubnames p,
         lateral (select name
                     from cities c
                order by c.pos <-> p.pos
                   limit 1) c
order by pos <-> point(-0.12,51.516) desc
   limit 5;

看來給定的 bounding box([bbox=-10.5,49.78,1.78,59])連愛爾蘭也涵蓋了……更重要的是,這個查詢的執行代價相當可觀:規劃器只會用「對每一筆 pubnames(27,878 次迴圈)都掃一次 cities」的方式解這個查詢,總共花了約 1.3 秒。

延伸輸出:慢查詢的 explain (analyze) 計畫
                                QUERY PLAN
════════════════════════════════════════════════════════════════════════
 Limit (actual time=1323.517..1323.518 rows=5 loops=1)
   -> Sort (actual time=1323.515..1323.515 rows=5 loops=1)
          Sort Key: ((p.pos <-> '(-0.12,51.516)'::point)) DESC
          Sort Method: top-N heapsort Memory: 25kB
          -> Nested Loop (actual time=0.116..1310.214 rows=27878 loops=1)
                -> Seq Scan on pubnames p (actual time=0.015..4.465 rows=27878 loops=1)
                -> Limit (actual time=0.044..0.044 rows=1 loops=27878)
                      -> Sort (actual time=0.043..0.043 rows=1 loops=27878)
                            Sort Key: ((c.pos <-> p.pos))
                            -> Seq Scan on cities c (... rows=73 loops=27878)
 Planning time: 0.236 ms
 Execution time: 1323.592 ms
(13 rows)

可以用 CTE 強迫規劃器照「顯而易見」的順序做:先取出最遠的五間酒吧,對這五間找最近城市。結果相同,但只花約 60ms,而不是超過一秒。

with pubs as (
    select name, pos,
           round((pos <@> point(-0.12,51.516))::numeric, 3) as miles
      from pubnames
  order by pos <-> point(-0.12,51.516) desc
     limit 5
)
select c.name as city, p.name, p.miles
  from pubs p, lateral (select name
                          from cities c
                      order by c.pos <-> p.pos
                         limit 1) c;

各城市最常見的酒吧名#

先找酒吧數最多的城市。由於 OSM 匯出的資料裡,一個城市只有一個代表點,我們的歸屬規則是:酒吧距城市位置點 5 英里以內就算屬於該城市——等於以該點畫一個 5 英里的圓,圓內的都算城裡的:

  select c.name, count(cp)
    from cities c, lateral (select name
                              from pubnames p
                              where (p.pos <@> c.pos) < 5) as cp
group by c.name
order by count(cp) desc
   limit 10;
    name     │ count
═════════════╪═══════
 London      │ 1388
 Westminster │ 1383
 Dublin      │   402
 Manchester  │   306
 Bristol     │   292
 Leeds       │   292
 Edinburgh   │   286
 Liverpool   │   258
 Nottingham  │   218
 Glasgow     │   217
(10 rows)

看地圖便知,依我們「5 英里內」的任意規則,Westminster 其實落在 London 之內,所以下個查詢直接把它濾掉。

留給讀者的練習:寫一個查詢,把實際上在 Westminster(距其位置 1 英里內)的酒吧從 London 的計數中扣除,再推廣到整個資料集的類似情況。另外提示:若應用需要考慮城市的真實形狀而非像這裡用猜的,該用 PostGIS 了。

最後是各城市最常見的酒吧名。計數時再次正規化酒吧名稱,但顯示時保留所有被計入的拼法:

  select c.name,
         array_to_string(array_agg(distinct(cp.name) order by cp.name), ', '),
         count(*)
    from cities c,
         lateral (select name
                    from pubnames p
                   where (p.pos <@> c.pos) < 5) as cp
   where c.name <> 'Westminster'
group by c.name, replace(replace(cp.name, 'The ', ''), 'And', '&')
order by count(*) desc
   limit 10;

這個查詢用上了前面所有的技巧:

  • lateral 子查詢
  • 在查詢內完成資料正規化
  • 用 earthdistance 提供的 <@> point 運算子計算距離
  • 以有序聚合去除重複值
   name   │            array_to_string             │ count
══════════╪════════════════════════════════════════╪═══════
 London   │ Prince of Wales, The Prince of Wales   │   15
 London   │ All Bar One                            │   12
 London   │ The Beehive                            │    8
 London   │ O'Neills                               │    7
 London   │ The Crown                              │    7
 London   │ The Windmill                           │    7
 London   │ Coach and Horses, The Coach and Horses │    6
 London   │ The Ship                               │    6
 Bradford │ New Inn, The New Inn                   │    6
 London   │ Red Lion, The Red Lion                 │    6
(10 rows)