PostGIS 擴充了 PostgreSQL 的功能，使其具備處理空間資料的能力，對於地理資訊系統和位置服務應用至關重要。本文的範例涵蓋了從農夫市場資料函式庫建立地理欄位、利用空間索引提升查詢效能，到計算市場間距離、尋找最近市場等空間分析技巧。此外，文章也詳細說明瞭如何匯入和處理 Shapefile 格式的美國人口普查資料，包含使用 shp2pgsql 工具和圖形介面，以及如何利用空間函式計算行政區面積等實務操作。這些技術為開發者提供了處理和分析空間資料的有效工具，可應用於各種地理資訊相關的應用場景。
使用PostGIS進行空間資料分析
PostGIS是一種功能強大的空間資料函式庫擴充套件，能夠讓我們在PostgreSQL資料函式庫中進行空間資料的儲存、查詢和分析。在本章中，我們將探討如何使用PostGIS來分析農夫市場的資料。
建立和填入地理欄位
要進行空間查詢，我們需要將經度和緯度座標轉換成一個具有空間資料型別的欄位。由於我們正在處理跨越整個美國的資料，因此我們需要使用geography型別來確保在大範圍內的距離測量是準確的。
ALTER TABLE farmers_markets ADD COLUMN geog_point geography(POINT,4326);
UPDATE farmers_markets
SET geog_point = ST_SetSRID(ST_MakePoint(longitude,latitude),4326)::geography;
CREATE INDEX market_pts_idx ON farmers_markets USING GIST (geog_point);
內容解密：
ALTER TABLE陳述式用於新增一個名為geog_point的欄位，其資料型別為geography(POINT,4326)，表示這是一個參考WGS 84座標系統（SRID 4326）的點。UPDATE陳述式用於填入geog_point欄位。ST_MakePoint(longitude,latitude)函式根據經度和緯度欄位建立一個幾何點，然後使用ST_SetSRID函式設定其SRID為4326，最後將其轉換為geography型別以符合geog_point欄位的型別。CREATE INDEX陳述式用於在geog_point欄位上建立一個GIST索引，以加速空間查詢。
分析農夫市場資料
美國農業部（USDA）的全國農夫市場目錄（National Farmers’ Market Directory）收錄了超過8,600個農夫市場的資訊，包括它們的位置和提供的產品。我們可以使用SQL空間查詢來找出最近的農夫市場。
首先，我們需要建立和填入一個包含農夫市場資料的表格。
CREATE TABLE farmers_markets (
fmid bigint PRIMARY KEY,
market_name text NOT NULL,
street text,
city text,
county text,
st text NOT NULL,
zip text,
longitude numeric(10,7),
latitude numeric(10,7),
organic text NOT NULL
);
COPY farmers_markets
FROM 'C:\YourDirectory\farmers_markets.csv'
WITH (FORMAT CSV, HEADER);
內容解密：
CREATE TABLE陳述式定義了farmers_markets表格的結構，包括欄位名稱和資料型別。COPY陳述式用於從CSV檔案中匯入資料到farmers_markets表格中。
空間索引的重要性
在進行空間查詢之前，建立空間索引是非常重要的。空間索引可以大幅提高查詢效能，尤其是在處理大量資料時。
CREATE INDEX market_pts_idx ON farmers_markets USING GIST (geog_point);
內容解密：
CREATE INDEX陳述式用於在geog_point欄位上建立一個GIST索引。- GIST索引是一種適合於空間資料的索引型別，可以有效地加速空間查詢。
檢視地理資料
建立好空間索引後，我們可以使用SELECT陳述式來檢視地理資料。
SELECT longitude, latitude, geog_point, ST_AsEWKT(geog_point)
FROM farmers_markets
WHERE longitude IS NOT NULL
LIMIT 5;
內容解密：
SELECT陳述式用於選擇特定的欄位，包括經度、緯度、geog_point和其擴充套件WKT表示。ST_AsEWKT(geog_point)函式用於將geog_point轉換為擴充套件WKT格式，以便於檢視其座標和SRID。
現在，我們已經準備好對這些點進行計算。接下來，我們將探討如何使用PostGIS進行更複雜的空間分析。
在給定距離內尋找地理位置
幾年前，在報導愛荷華州農業相關的新聞時，我曾造訪德梅因市（Des Moines）的Downtown Farmers’ Market。這座大型農夫市場匯集了數百家攤販，橫跨了愛荷華州首府的數個市區街區。農業在當地是一項重要的產業，即使是市中心的農夫市場也相當龐大，但它並不是該地區唯一的農夫市場。讓我們利用PostGIS來找出更多鄰近市中心的其他農夫市場。
使用ST_DWithin()函式
PostGIS的ST_DWithin()函式能夠傳回一個布林值，指示一個空間物件是否在另一個物件的指定距離內。如果使用的是地理（geography）資料型別，就必須以公尺作為距離單位；如果是幾何（geometry）型別，則使用由SRID指定的距離單位。
注意：PostGIS的距離測量對於幾何資料是直線距離，而對於地理資料則是在球面上計算。請注意不要將其與道路駕駛距離混淆，後者通常會更遠。要進行與駕駛距離相關的計算，可以參考pgRouting擴充套件套件，網址為https://pgrouting.org/。
範例程式碼
清單15-10使用了ST_DWithin()函式來篩選農夫市場，顯示出在德梅因市Downtown Farmers’ Market周圍10公里內的市場。
SELECT market_name,
city,
st,
geog_point
FROM farmers_markets
WHERE ST_DWithin(geog_point,
ST_GeogFromText('POINT(-93.6204386 41.5853202)'),
10000)
ORDER BY market_name;
內容解密：
geog_point是 farmers_markets 表中儲存市場位置的地理資料欄位。ST_GeogFromText('POINT(-93.6204386 41.5853202)')將WKT格式的座標轉換為地理資料型別的點，代表Downtown Farmers’ Market的位置。10000是距離單位（公尺），等於10公里。
執行結果
執行查詢後，應傳回9行結果（省略了geog_point欄位以節省篇幅）：
market_name city st
---
-
---
-
---
-
---
-
---
-
---
-
---
-
---
-
---
--
---
-
---
-
---
---
-
---
-
Beaverdale Farmers Market Des Moines Iowa
Capitol Hill Farmers Market Des Moines Iowa
Downtown Farmers' Market - Des Moines Des Moines Iowa
Drake Neighborhood Farmers Market Des Moines Iowa
Eastside Farmers Market Des Moines Iowa
Highland Park Farmers Market Des Moines Iowa
Historic Valley Junction Farmers Market West Des Moines Iowa
LSI Global Greens Farmers' Market Des Moines Iowa
Valley Junction Farmers Market West Des Moines Iowa
這些結果顯示了鄰近德梅因市Downtown Farmers’ Market的其他農夫市場。
計算地理位置之間的距離
要知道這些市場與市中心之間的確切距離，可以使用ST_Distance()函式。
使用ST_Distance()函式
ST_Distance()函式傳回兩個幾何物件之間的最小距離，對於地理資料傳回公尺，對於幾何資料傳回SRID單位。
範例程式碼
清單15-11計算了洋基體育場（Yankee Stadium）與梅斯棒球場（Citi Field）之間的距離（英里）。
SELECT ST_Distance(
ST_GeogFromText('POINT(-73.9283685 40.8296466)'),
ST_GeogFromText('POINT(-73.8480153 40.7570917)')
) / 1609.344 AS mets_to_yanks;
內容解密：
ST_GeogFromText()將WKT格式的座標轉換為地理資料型別的點，分別代表洋基體育場和梅斯棒球場的位置。ST_Distance()計算兩個點之間的距離，傳回公尺。- 除以
1609.344將公尺轉換為英里。
應用ST_Distance()於農夫市場資料
清單15-12再次找出德梅因市Downtown Farmers’ Market周圍10公里內的農夫市場，並顯示距離（英里）。
SELECT market_name,
city,
round(
(ST_Distance(geog_point,
ST_GeogFromText('POINT(-93.6204386 41.5853202)')) / 1609.344)::numeric, 2
) AS miles_from_dt
FROM farmers_markets
WHERE ST_DWithin(geog_point,
ST_GeogFromText('POINT(-93.6204386 41.5853202)'),
10000)
ORDER BY miles_from_dt ASC;
內容解密：
round()函式用於將距離四捨五入到小數點後兩位。ST_Distance()計算每個市場與Downtown Farmers’ Market之間的距離。WHERE子句使用ST_DWithin()篩選出在10公里內的市場。
結果展示
查詢結果按照距離升序排列：
market_name city miles_from_dt
---
-
---
-
---
-
---
-
---
-
---
-
---
-
---
-
---
--
---
-
---
-
---
---
-
---
-
---
-
---
--
Downtown Farmers' Market - Des Moines Des Moines 0.00
Capitol Hill Farmers Market Des Moines 1.15
Drake Neighborhood Farmers Market Des Moines 1.70
LSI Global Greens Farmers' Market Des Moines 2.30
Highland Park Farmers Market Des Moines 2.93
Eastside Farmers Market Des Moines 3.40
Beaverdale Farmers Market Des Moines 3.74
Historic Valley Junction Farmers Market West Des Moines 4.68
透過這些範例，我們展示瞭如何使用PostGIS的ST_DWithin()和ST_Distance()函式來找出鄰近的地理位置並計算它們之間的距離。這些技術對於需要進行空間分析和查詢的應用程式非常有用。
使用PostGIS進行地理空間分析
尋找最近的地理物件
有時，我們需要資料函式庫傳回最接近另一個物件的空間物件，而無需指定搜尋的任意距離。例如，我們可能想要找到最近的農貿市場，無論它距離10公里還是100公里。要做到這一點，我們可以指示PostGIS使用查詢的ORDER BY子句中的<->距離運算子實作K-最近鄰搜尋演算法。
K-最近鄰搜尋示例
假設我們計劃存取緬因州巴爾港的度假勝地，並希望找到離城鎮最近的三個農貿市場。我們可以使用清單15-13中的程式碼。
SELECT market_name,
city,
st,
round(
(ST_Distance(geog_point,
ST_GeogFromText('POINT(-68.2041607 44.3876414)')
) / 1609.344)::numeric, 2
) AS miles_from_bh
FROM farmers_markets
ORDER BY geog_point <-> ST_GeogFromText('POINT(-68.2041607 44.3876414)')
LIMIT 3;
內容解密：
SELECT market_name, city, st：選擇農貿市場的名稱、城市和州。ST_Distance(geog_point, ST_GeogFromText('POINT(-68.2041607 44.3876414)'))：計算每個農貿市場與巴爾港市中心的距離。round(...::numeric, 2)：將距離四捨五入到小數點後兩位，並轉換為英里。ORDER BY geog_point <-> ST_GeogFromText('POINT(-68.2041607 44.3876414)')：按距離排序結果，使用<->運算子實作K-最近鄰搜尋。LIMIT 3：限制結果為最近的三個農貿市場。
處理人口普查Shapefiles
Shapefile是一種由Esri開發的GIS資料檔案格式，用於儲存具有地理特徵的資料。Shapefiles包含描述特徵形狀的資訊，以及每個特徵的屬性資料函式庫。這些屬性可能包括名稱和其他人口統計描述符。
Shapefile的結構和內容
一個Shapefile由多個具有不同副檔名的檔案組成，每個檔案都有不同的用途。常見的副檔名包括：
.shp：儲存特徵幾何的主檔案。.shx：儲存特徵幾何索引的索引檔案。.dbf：儲存特徵屬性資訊的資料函式庫表（以dBASE格式）。.xml：儲存Shapefile元資料的XML格式檔案。.prj：儲存坐標系資訊的投影檔案。
載入Shapefiles
如果您使用的是Windows，PostGIS套件包含一個具有簡單圖形使用者介面（GUI）的Shapefile匯入/匯出管理器。在macOS和Linux發行版上，我們將使用命令列應用程式shp2pgsql。
Windows Shapefile匯入器/匯出器
在Windows上，如果您按照第1章中的安裝步驟進行操作，則應該可以透過選擇“開始”▶“PostGIS Bundle x.y for PostgreSQL x64 x.y”▶“PostGIS Bundle x.y for PostgreSQL x64 x.y Shapefile and DBF Loader Exporter”找到Shapefile匯入/匯出管理器。點選以啟動應用程式。
探索和分析Shapefiles
一旦載入Shapefile，您就可以存取其空間物件和每個物件的屬性。我們將在後面的章節中探討一些額外的分析功能。
將應用程式與分析資料函式庫連線
要建立應用程式與分析資料函式庫之間的連線，請依照以下步驟進行：
設定PostGIS連線
- 點選「檢視連線詳細資訊」。
- 在開啟的對話方塊中，輸入
postgres作為使用者名稱，如果在初始設定時為伺服器新增了密碼，請輸入該密碼。 - 確認伺服器主機（Server Host）預設為
localhost，連線埠（Port）預設為5432。除非連線到不同的伺服器或連線埠，否則保持預設值不變。 - 輸入
analysis作為資料函式庫名稱。圖15-4展示了連線設定的範例。
載入Shapefile
- 點選「確定」。在日誌視窗中應看到「連線成功」的訊息。成功建立PostGIS連線後，即可載入Shapefile。
- 在「選項」下，將DBF檔案字元編碼更改為
Latin1，因為Shapefile屬性包含需要此編碼的縣名稱。保持預設的核取方塊，包括在空間欄位上建立索引的選項。點選「確定」。 - 點選「新增檔案」並選擇之前儲存的
tl_2019_us_county.shp。點選「開啟」。該檔案應出現在載入器的Shapefile清單中，如圖15-5所示。 - 在「表格」欄中，雙擊以選擇表格名稱，將其替換為
us_counties_2019_shp。按下ENTER鍵接受該值。 - 在「SRID」欄中，雙擊並輸入
4269。這是美國聯邦機構（包括美國人口普查局）常用的北美1983座標系統的ID。再次按下ENTER鍵接受該值。 - 點選「匯入」。
驗證匯入結果
在日誌視窗中，應看到一則訊息，最後顯示：
Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Shapefile import completed.
切換到pgAdmin，在物件瀏覽器中展開analysis節點，並繼續展開，選擇Schemas▶public▶Tables。右鍵點選Tables並選擇「重新整理」，應看到us_counties_2019_shp列出。恭喜！已成功將Shapefile載入表格。
使用shp2pgsql匯入Shapefile
對於macOS和Linux上的某些PostGIS發行版，Shapefile匯入/匯出管理器不可用。因此，將展示如何使用PostGIS命令列工具shp2pgsql匯入Shapefile。
命令列語法
shp2pgsql -I -s SRID -W encoding shapefile_name table_name | psql -d database -U user
-I在新表格的幾何欄位上使用GiST新增索引。-s指定幾何資料的SRID。-W指設定檔案編碼（如有必要）。shapefile_name是副檔名為.shp的檔案名稱（包含完整路徑）。table_name是要匯入Shapefile的新表格名稱。
範例命令
shp2pgsql -I -s 4269 -W LATIN1 tl_2019_us_county.shp us_counties_2019_shp | psql -d analysis -U postgres
探索Census 2019 Counties Shapefile
新的us_counties_2019_shp表格包含各縣的名稱以及聯邦資訊處理標準（FIPS）程式碼等欄位。geom欄位包含每個縣界的空間資料。
檢查geom欄位的WKT表示法
SELECT ST_AsText(geom)
FROM us_counties_2019_shp
ORDER BY gid
LIMIT 1;
結果顯示一個具有數百個座標對的MultiPolygon。
找出面積最大的縣
使用ST_Area()函式計算面積
SELECT name,
statefp AS st,
round(
(ST_Area(geom::geography) / 2589988.110336)::numeric, 2
) AS square_miles
FROM us_counties_2019_shp
ORDER BY square_miles DESC
LIMIT 5;
此查詢計算每個縣的面積（單位：平方英里），並找出面積最大的前5個縣。注意，ST_Area()函式在幾何資料型別上傳回SRID單位的值，因此需要將其轉換為地理資料型別以獲得以平方米為單位的值，然後再轉換為平方英里。
