GeoNames 地理資料庫涵蓋所有國家,收錄超過一千一百萬筆地名,可免費下載。
GeoNames 網站提供線上查詢,所有資料也開放下載使用。一如往常,資料是 ad-hoc 格式,得先經過處理與正規化才能在 PostgreSQL 裡好好使用。
做法是先建立 raw schema,把官方發佈的檔案原樣載入:raw.geonames(主要地名資料,含 geonameid、name、經緯度、feature_class/feature_code、country_code、admin1–4 code、population 等欄位)、raw.country、raw.feature、raw.admin1、raw.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是最佳選擇(唯一且為整數)。 languages與neighbours都是多值欄位:逗號分隔的語言或國碼清單。- 要達到 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_code、admin1_code、admin2_code 做地理分層。以法國資料為例,admin1_code/admin2_code 是 44、67 之類的代碼,必須展開才有意義(美國資料的 admin1 是 IL 這種州代碼,正規化的必要性比較不明顯)。
探索資料時可以用
offset隨手翻頁,但如前面章節所述,應用程式查詢絕對不要用 offset。這裡是互動式資料探索,才勉強可接受。
GeoNames 提供 admin1CodesASCII.txt 與 admin2Codes.txt 供正規化使用,代碼同樣是 AD.06、AF.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 列——畢竟取樣是隨機演算法。