GeoNames 地理資料庫涵蓋所有國家,收錄超過一千一百萬筆地名,可免費下載。

GeoNames 網站提供線上查詢,所有資料也開放下載使用。一如往常,資料是 ad-hoc 格式,得先經過處理與正規化才能在 PostgreSQL 裡好好使用。

做法是先建立 raw schema,把官方發佈的檔案原樣載入:raw.geonames(主要地名資料,含 geonameid、name、經緯度、feature_class/feature_code、country_code、admin1–4 code、population 等欄位)、raw.countryraw.featureraw.admin1raw.admin2,多數檔案用 \copy ... with csv delimiter E'\t' 即可載入。

raw schema 載入腳本(節錄)
begin;

create schema if not exists raw;

create table raw.geonames
 (
   geonameid         bigint,
   name              text,
   asciiname         text,
   alternatenames    text,
   latitude          double precision,
   longitude         double precision,
   feature_class     text,
   feature_code      text,
   country_code      text,
   cc2               text,
   admin1_code       text,
   admin2_code       text,
   admin3_code       text,
   admin4_code       text,
   population        bigint,
   elevation         bigint,
   dem               bigint,
   timezone          text,
   modification      date
 );

create table raw.country
 (
   iso                text,
   iso3               text,
   isocode            integer,
   fips               text,
   name               text,
   capital            text,
   area               double precision,
   population         bigint,
   continent          text,
   tld                text,
   currency_code      text,
   currency_name      text,
   phone              text,
   postal_code_format text,
   postal_code_regex  text,
   languages          text,
   geonameid          bigint,
   neighbours         text,
   fips_equiv         text
 );

\copy raw.country from 'countryInfoData.txt' with csv delimiter E'\t'

create table raw.feature
 (
   code        text,
   description text,
   comment     text
 );

\copy raw.feature from 'featureCodes_en.txt' with csv delimiter E'\t'

create table raw.admin1
 (
   code       text,
   name       text,
   ascii_name text,
   geonameid  bigint
 );

\copy raw.admin1 from 'admin1CodesASCII.txt' with csv delimiter E'\t'

create table raw.admin2
 (
   code       text,
   name       text,
   ascii_name text,
   geonameid  bigint
 );

\copy raw.admin2 from 'admin2Codes.txt' with csv delimiter E'\t'

commit;

唯獨 raw.geonames 沒有 \copy 指令:因為部分地名含未正確引用、也未跳脫的單雙引號,copy 會失敗。改用 pgloader 載入:

load csv
  from /tmp/geonames/allCountries.txt
  into pgsql://appdev@/appdev
  target table raw.geonames

  with fields terminated by '\t',
       fields optionally enclosed by '§',
       fields escaped by '%',
       truncate;

在準備本書的筆電上,11,540,466 列、1.5 GB 的資料約 6 分 43 秒載入完成。

原始資料載好後即可正規化。基本原則:避免屬性之間的依賴——一旦有依賴,就把那組自成一體的資料拆到獨立的表管理。raw.geonames 用到多份 GeoNames 另外提供的參照資料,所以先從參照資料下手。

特徵分類#

GeoNames 為所有地理資料標上 feature class 與 feature 兩層代碼,說明在 featureCodes_en.txt;class 的文字說明只有網頁上有,得手動整理。

begin;

create schema if not exists geoname;

create table geoname.class
 (
   class       char(1) not null primary key,
   description text
 );

insert into geoname.class (class, description)
     values ('A', 'country, state, region,...'),
            ('H', 'stream, lake, ...'),
            ('L', 'parks,area, ...'),
            ('P', 'city, village,...'),
            ('R', 'road, railroad '),
            ('S', 'spot, building, farm'),
            ('T', 'mountain,hill,rock,... '),
            ('U', 'undersea'),
            ('V', 'forest,heath,...');

create table geoname.feature
 (
   class       char(1) not null references geoname.class(class),
   feature     text    not null,
   description text,
   comment     text,

  primary key(class, feature)
 );

insert into geoname.feature
     select substring(code from 1 for 1) as class,
            substring(code from 3) as feature,
            description,
            comment
       from raw.feature
      where feature.code <> 'null';

commit;
  • 檔案末行有一筆字面上就是四個字母 null 的項目,必須排除。
  • 原始檔用 A.ADM1 這種寫法表示 class A、feature ADM1;正規化時拆成兩個屬性,自然鍵就是 class + feature 的組合。
資料載好後的 top-10 統計
  select class, feature, description, count(*)
    from feature
         left join geoname using(class,feature)
group by class, feature
order by count desc
   limit 10;
 class │ feature │   description   │  count
═══════╪═════════╪═════════════════╪═════════
 P     │ PPL     │ populated place │ 1711458
 H     │ STM     │ stream          │  300283
 S     │ CH      │ church          │  236394
 S     │ FRM     │ farm            │  234536
 S     │ SCH     │ school          │  223402
 T     │ HLL     │ hill            │  212659
 T     │ MT      │ mountain        │  192454
 S     │ HTL     │ hotel           │  170896
 H     │ LK      │ lake            │  162922
 S     │ BLDG    │ building(s)     │  143742
(10 rows)

國家#

raw.country(每列含 iso、iso3、isocode、name、capital、continent、tld、currency、languages、neighbours 等欄位)有幾個明顯的正規化問題:

  • 沒有任何機制保證無重複列,需要加主鍵——isocode 是最佳選擇(唯一且為整數)。
  • languagesneighbours 都是多值欄位:逗號分隔的語言或國碼清單。
  • 要達到 2NF,所有非鍵屬性必須依賴整個鍵,但貨幣與郵遞區號格式並不依賴於國家。

檢查屬性是否真的依賴鍵,可以用這類查詢:

  select currency_code, currency_name, count(*)
    from raw.country
group by currency_code, currency_name
order by count desc
   limit 5;

結果顯示有 34 個國家共用 Euro——貨幣顯然不是國家的函數。

本書略過貨幣、語言與郵遞區號格式,只保留部分資訊。正規化過程:

begin;

create schema if not exists geoname;

create table geoname.continent
 (
   code   char(2) primary key,
   name   text
 );

insert into geoname.continent(code, name)
     values ('AF', 'Africa'),
            ('NA', 'North America'),
            ('OC', 'Oceania'),
            ('AN', 'Antarctica'),
            ('AS', 'Asia'),
            ('EU', 'Europe'),
            ('SA', 'South America');

create table geoname.country
 (
   isocode   integer primary key,
   iso       char(2) not null,
   iso3      char(3) not null,
   fips      text,
   name      text,
   capital   text,
   continent char(2) references geoname.continent(code),
   tld       text,
   geonameid bigint
 );

insert into geoname.country
     select isocode, iso, iso3, fips, name,
            capital, continent, tld, geonameid
       from raw.country;

create table geoname.neighbour
 (
   isocode   integer not null references geoname.country(isocode),
   neighbour integer not null references geoname.country(isocode),

  primary key(isocode, neighbour)
 );

insert into geoname.neighbour
   with n as(
     select isocode,
            regexp_split_to_table(neighbours, ',') as neighbour
       from raw.country
   )
   select n.isocode,
          country.isocode
     from n
          join geoname.country
            on country.iso = n.neighbour;

commit;

補上洲別清單(讓地區下鑽更完整),並把多值的 neighbours 拆成 geoname.neighbour 關聯表——每個國家與其接壤鄰國的關係,之後查詢就很直接:

select neighbour.iso,
       neighbour.name,
       neighbour.capital,
       neighbour.tld

  from geoname.neighbour as border

          join geoname.country as country
            on border.isocode = country.isocode

          join geoname.country as neighbour
            on border.neighbour = neighbour.isocode

 where country.iso = 'FR';
 iso │    name     │      capital     │ tld
═════╪═════════════╪══════════════════╪═════
 CH  │ Switzerland │ Bern             │ .ch
 DE  │ Germany     │ Berlin           │ .de
 BE  │ Belgium     │ Brussels         │ .be
 LU  │ Luxembourg  │ Luxembourg       │ .lu
 IT  │ Italy       │ Rome             │ .it
 AD  │ Andorra     │ Andorra la Vella │ .ad
 MC  │ Monaco      │ Monaco           │ .mc
 ES  │ Spain       │ Madrid           │ .es
(8 rows)

行政區劃#

原始資料以 country_codeadmin1_codeadmin2_code 做地理分層。以法國資料為例,admin1_codeadmin2_code4467 之類的代碼,必須展開才有意義(美國資料的 admin1 是 IL 這種州代碼,正規化的必要性比較不明顯)。

探索資料時可以用 offset 隨手翻頁,但如前面章節所述,應用程式查詢絕對不要用 offset。這裡是互動式資料探索,才勉強可接受。

GeoNames 提供 admin1CodesASCII.txtadmin2Codes.txt 供正規化使用,代碼同樣是 AD.06AF.01.1125426 這種點分格式,正好在此拆開並加上約束以確保資料品質:

begin;

create schema if not exists geoname;

create table geoname.region
 (
   isocode   integer not null references geoname.country(isocode),
   regcode   text not null,
   name      text,
   geonameid bigint,

  primary key(isocode, regcode)
 );

insert into geoname.region
   with admin as
   (
     select regexp_split_to_array(code, '[.]') as code,
            name,
            geonameid
       from raw.admin1
   )
   select country.isocode as isocode,
          code[2] as regcode,
          admin.name,
          admin.geonameid
     from admin
          join geoname.country
            on country.iso = code[1];

create table geoname.district
 (
   isocode   integer not null,
   regcode   text not null,
   discode   text not null,
   name      text,
   geonameid bigint,

  primary key(isocode, regcode, discode),
  foreign key(isocode, regcode)
    references geoname.region(isocode, regcode)
 );

insert into geoname.district
   with admin as
   (
     select regexp_split_to_array(code, '[.]') as code,
            name,
            geonameid
       from raw.admin2
   )
     select region.isocode,
            region.regcode,
            code[3],
            admin.name,
            admin.geonameid
       from admin
            join geoname.country
              on country.iso = code[1]
            join geoname.region
              on region.isocode = country.isocode
             and region.regcode = code[2];

commit;

之前的查詢現在可以改寫成顯示地區與行政區名稱(內部仍保有代碼備用)。查詢用 left join,因為有些地理資料缺 admin1 或 admin2 層級的細節:

        name         │  region   │           district
═════════════════════╪═══════════╪═══════════════════════════════
 Zintzel du Nord     │ Grand Est │ Département du Bas-Rhin
 Zinswiller          │ Grand Est │ Département du Bas-Rhin
 Ruisseau de Zingajo │ Corsica   │ Département de la Haute-Corse
 Zincourt            │ Grand Est │ Département des Vosges
 Zimming             │ Grand Est │ Département de la Moselle
(5 rows)

地理定位資料#

載入主資料前,先研究資料的分層情形:

select count(*) as all,
       count(*) filter(where country_code is null) as no_country,
       count(*) filter(where admin1_code is null) as no_region,
       count(*) filter(where admin2_code is null) as no_district,
       count(*) filter(where feature_class is null) as no_class,
       count(*) filter(where feature_code is null) as no_feat
  from raw.geonames;
   all    │ no_country │ no_region │ no_district │ no_class │ no_feat
══════════╪════════════╪═══════════╪═════════════╪══════════╪═════════
 11540466 │       5821 │     45819 │     5528455 │     5074 │   95368
(1 row)

不少項目沒有國家參照,更多項目缺細部區劃(admin1/admin2 不一定存在),也有無 feature 與 class 的點,有些位於北極圈。

因此正規化查詢必須使用 left join,允許外鍵參照不存在時欄位為 null;並且要按「國家 → 地區 → 行政區」的順序逐層下鑽,因為資料集包含多種精度層級的點。

begin;

create table geoname.geoname
 (
   geonameid         bigint primary key,
   name              text,
   location          point,
   isocode           integer,
   regcode           text,
   discode           text,
   class             char(1),
   feature           text,
   population        bigint,
   elevation         bigint,
   timezone          text,

   foreign key(isocode)
    references geoname.country(isocode),

   foreign key(isocode, regcode)
    references geoname.region(isocode, regcode),

   foreign key(isocode, regcode, discode)
    references geoname.district(isocode, regcode, discode),

   foreign key(class)
    references geoname.class(class),

   foreign key(class, feature)
    references geoname.feature(class, feature)
 );

insert into geoname.geoname
  with geo as
  (
      select geonameid,
             name,
             point(longitude, latitude) as location,
             country_code,
             admin1_code,
             admin2_code,
             feature_class,
             feature_code,
             population,
             elevation,
             timezone
        from raw.geonames
    )
      select geo.geonameid,
             geo.name,
             geo.location,
             country.isocode,
             region.regcode,
             district.discode,
             feature.class,
             feature.feature,
             population,
             elevation,
             timezone
        from geo
             left join geoname.country
               on country.iso = geo.country_code

             left join geoname.region
               on region.isocode = country.isocode
              and region.regcode = geo.admin1_code

             left join geoname.district
               on district.isocode = country.isocode
              and district.regcode = geo.admin1_code
              and district.discode = geo.admin2_code

             left join geoname.feature
               on feature.class = geo.feature_class
              and feature.feature = geo.feature_code;

create index on geoname.geoname using gist(location);

commit;

注意經緯度合併成了 PostgreSQL 的 point 型別,並建立 GiST 索引。載好後即可輕鬆分析各洲分佈——資料明顯偏向亞洲、北美與歐洲,南極洲則相當稀疏:

     name      │  count  │  pct  │               hist
═══════════════╪═════════╪═══════╪═══════════════════════════════════
 Africa        │ 1170043 │ 10.14 │ ■■■■■■■■■■
 Antarctica    │   21125 │  0.18 │
 Asia          │ 3772195 │ 32.70 │ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
 Europe        │ 2488807 │ 21.58 │ ■■■■■■■■■■■■■■■■■■■■■■
 North America │ 3210802 │ 27.84 │ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■
 Oceania       │  354325 │  3.07 │ ■■■
 South America │  517347 │  4.49 │ ■■■■
(7 rows)

地理定位的 GiST 索引#

gist(location) 索引在 1,150 萬筆資料中搜尋特定位置時非常有用。PostgreSQL 支援多種索引掃描情境,包括 kNN 查找(最近鄰查找,nearest neighbor lookup)。

先前在陣列型別範例中載入過 20 萬筆帶地理位置的推文(hashtag 表)。透過 lateral left join 加上 order by ... <-> ... limit k,就能為每則推文找出最近的 GeoNames 地點,再接上正規化好的國家、地區、行政區資訊:

  select id,
         round((hashtag.location <-> geoname.location)::numeric, 3) as dist,
         country.iso,
         region.name as region,
         district.name as district
     from hashtag
         left join lateral
         (
              select geonameid, isocode, regcode, discode, location
                from geoname.geoname
            order by location <-> hashtag.location
               limit 1
         )
         as geoname
         on true
         left join geoname.country using(isocode)
         left join geoname.region using(isocode, regcode)
         left join geoname.district using(isocode, regcode, discode)
order by id
   limit 5;

<-> 運算子計算兩點距離,配合 limit 1 便為每筆 hashtag 選出最近的已知地點:

         id         │ dist  │ iso │    region    │      district
════════════════════╪═══════╪═════╪══════════════╪═════════════════════
 720553447402160128 │ 0.004 │ US  │ Florida      │ Orange County
 720553457015324672 │ 0.004 │ US  │ Texas        │ Smith County
 720553458596757504 │ 0.001 │ US  │ Florida      │ Orange County
 720553466804989952 │ 0.001 │ US  │ Pennsylvania │ Philadelphia County
 720553475923271680 │ 0.000 │ US  │ New York     │ Nassau County
(5 rows)

explain (costs off) 檢查可看到 Index Scan using geoname_location_idx on geoname——GiST 索引確實被使用,整個查詢在筆電上約 13 毫秒完成。

完整查詢計畫(explain costs off)
                              QUERY PLAN
══════════════════════════════════════════════════════════════════════
 Limit
   -> Nested Loop Left Join
        -> Nested Loop Left Join
             -> Nested Loop Left Join
                   Join Filter: (geoname.isocode = country.isocode)
                   -> Nested Loop Left Join
                         -> Index Scan using hashtag_pkey on hashtag
                         -> Limit
                               -> Index Scan using geoname_location_idx on geoname
                                     Order By: (location <-> hashtag.location)
                   -> Materialize
                         -> Seq Scan on country
             -> Index Scan using region_pkey on region
                   Index Cond: ((geoname.isocode = isocode) AND (geoname.regcode = regcode))
        -> Index Scan using district_pkey on district
             Index Cond: ((geoname.isocode = isocode) AND (geoname.regcode = regcode)
                          AND (geoname.discode = discode))
(16 rows)

國家取樣#

1,100 多萬列的資料集不便隨書附上(Full/Enterprise Edition 附資料庫 dump 或 Docker 映像),因此改附 1% 的隨機樣本,用 PostgreSQL 的 tablesample 功能製作:

begin;

create schema if not exists sample;

drop table if exists sample.geonames;

create table sample.geonames
   as select geonameid,
             name,
             longitude,
             latitude,
             feature_class,
             feature_code,
             country_code,
             admin1_code,
             admin2_code,
             population,
             elevation,
             timezone
        from raw.geonames TABLESAMPLE bernoulli(1);

\copy sample.geonames to 'allCountries.sample.copy'

commit;

PostgreSQL 內建兩種取樣方法,參數都是 0–100 的百分比:

  • BERNOULLI 掃描整張表,逐列以指定機率獨立決定取或不取。
  • SYSTEM 做區塊層級取樣,被選中的區塊整塊回傳;小比例取樣時明顯較快,但受叢集效應影響,樣本隨機性較差。

也可自行實作取樣方法,見文件 Writing A Table Sampling Method

執行結果約 115,904 列;再跑一次得到 115,071 列——畢竟取樣是隨機演算法。