我一直是一個地圖極客，可以追溯到 80 年代，當時我會拿一張道路地圖集和一些描圖紙，然後在我自己的道路網絡中繪制。 我最喜歡的遊戲之一是拿一張舊地圖或地球儀，并嘗試根據國家的名稱和形狀來确定它的制作時間。 在線上地圖軟體、大資料和資料科學時代，它變得更加有趣。 現在我可以下載下傳地圖的整個資料集并針對它編寫程式了！
推薦：用 NSDT設計器 快速搭建可程式設計3D場景。
1、Project Remote
最近，我一直在關注 Project Remote，這是一個家庭冒險計算和記錄美國每個州離道路最遠的點。 在撰寫本文時，他們已經完成了 33 個州，而俄勒岡州不是其中之一。 這讓我開始思考。 Project Remote 沒有公布他們的确切公式，但他們給出了一些提示。 基于這些，我可以計算出俄勒岡州最偏遠的點嗎？
好吧，這有點複雜。 俄勒岡州有成千上萬條道路，我如何計算那麼多數字？ 也許我會從小一點開始。
我将隻使用高速公路而不是所有道路，因為高速公路并不多。 我将從波特蘭開始，因為那是大多數高速公路所在的地方，也是人們熟悉的地方。
似乎有兩種方法可以找到距離高速公路最遠的點：
- 使用蠻力方法，設定點網格并計算每個點與最近的高速公路的距離
- 編寫一個 AI 從任意點開始并遠離高速公路
對于 Metal Toad 2017 年 12 月的資料科學黑客馬拉松，我和我的同僚 Kalina Wilson 認為 AI 選項會更有趣。
那麼，怎麼做呢？ 首先，我有以下問題：
- 給定一組緯度/經度坐标，計算機能否确定該點是否在俄勒岡州内？
- 計算機可以測量到最近的高速公路的距離嗎？
- 計算機能否将坐标移動到離高速公路較遠的附近點？
事實證明，所有這些問題的答案都是“是”，而且比我預期的要容易。
2、資料和工具
首先，我們需要一些資料和一些工具來分析它。 我選擇了以下工具，因為它們都是免費和開源的，而且它們都可以很好地協同工作。
2.1 地圖資料 - OpenStreetMap
OpenStreetMap 是開源地圖軟體。 其地圖資料完全由社群建構（即衆包），可免費下載下傳。 資料集包含州、縣、市的形狀資訊； 道路和高速公路； 江河湖海； 國家公園和森林； 以及企業、郵局和山峰等興趣點。
全球 OpenStreetMap 資料集壓縮為 40 GB。 這太啰嗦了，是以各種組織都提取了特定的區域，比如一個國家或州。 對于這個項目，我們使用了俄勒岡州的資料，即 134 MB。 你可以在 OpenStreetMap wiki 上閱讀所有不同的下載下傳選項。
2.2 GIS資料庫 - PostGIS
PostGIS 是 PostgreSQL 的擴充，用于處理幾何和地理資料。 它提供了一個 Geometry 列類型和幾十個與幾何相關的函數。 可用函數可以執行如下計算：
- 找出兩個形狀之間的距離
- 找到從一個形狀或點到另一個的方向
- 判斷一個點是否在多邊形内
- 機關轉換，例如，将緯度/經度轉換為米/英尺
OpenStreetMap 資料采用稱為 PBF 的特定格式，我們可以使用稱為 osm2pgsql 的開源工具将其加載到 PostgreSQL 中。
2.3 地理資訊系統 - QGIS
QGIS 是一個用于檢視地理資訊系統資料的桌面應用程式。它有一個可以直接連接配接到 PostgreSQL 的資料庫管理器。 你可以使用 PostGIS 函數和正常 SQL“where”子句編寫 SQL 查詢，并将結果作為圖層包含在地圖上。
2.4 地理編碼處理庫 - Geopy
Geopy 是一個 Python 庫，用于處理地理編碼資料并執行幾何和地理計算。 它可以通過調用各種網絡服務來查找位置和位址。 它還可以執行計算以測量兩點之間的距離。 使用一個幾乎沒有記錄的功能，我了解到它還可以做諸如“從（北緯 45.51°，西經 122.67°），以羅盤航向東北 45° 行駛 500 米”之類的事情。
2.5 前端可視化 - Leaflet|MapBox|Turf
為了在浏覽器中可視化，我們使用了：
- Leaflet JS ， 這是一個流行的 JavaScript 庫，用于建立互動式地圖。
- MapBox GL JS ，最初基于 Leaflet JS，但使用矢量切片，這增加了重要的用戶端功能（縮放、中繼資料查詢、樣式）。
- Turf JS ， 添加了空間分析和統計工具，可以進行測量、計算和選擇。
3、處理資料
将俄勒岡資料集加載到資料庫後，我可以使用資料庫前端工具檢視資料。 我發現我現在有了三個有趣的表：點、線和多邊形。 每個代表一種幾何形狀并包含具有該幾何形狀的每個對象。 “點”包含商戶、山峰等點源要素。 “線”包括道路、河流、湖泊和其他“開放”的幾何形狀。 “多邊形”表示閉合形狀，如州、縣、城市、國家森林和湖泊。
每個表都包含許多包含分類資訊（道路類型、水域類型、邊界類型、海拔等）的文本列和另一列稱為“way”的列。 “way”是 OpenStreetMap 所說的地理形狀。 它表示為 PostGIS“幾何”資料類型。 在資料庫中，Way 看起來像 0101000020E6100000FA7E6ABC74AB5EC0DF4F8D976E424740，它是地理形狀“POINT(-122.679 46.519)”的編碼表示。
PostGIS 提供了許多用于在 SQL 中操作方式和形狀的函數。
示例：查找從波特蘭先鋒法院廣場到時尚中心的方向和距離：
SELECT
ST_Distance(
ST_Transform(ST_GeomFromText('POINT(-122.679 45.519)', 4326), 2269),
ST_Transform(ST_GeomFromText('POINT(-122.667 45.532)', 4326), 2269)
) AS distance,
degrees(ST_Azimuth(
ST_GeomFromText('POINT(-122.679 45.519)', 4326),
ST_GeomFromText('POINT(-122.667 45.532)', 4326)
)) AS direction; 結果如下：
| Distance | Direction |
| 5650.219171193315 | 42.709389957366675 |
用人類的話說，這意味着 5650 英尺，航向為東北 42.7°。
再例如，我們可以使用函數 ST_Intersects() 來查找某個點所在的城市、縣和州的名稱。 請記住，市、縣和州邊界存儲在資料庫的“polygon”表中。 它們由“boundary”和“admin_level”字段辨別，其中“admin_level”為 4 是州，6 是縣，8 是市，10 是社群。
SELECT name
FROM planet_osm_polygon
WHERE
boundary = 'administrative'
AND ST_Intersects(
ST_GeomFromText('POINT(-122.679 45.519)', 4326),
ST_Transform(way, 4326)
) 結果如下：
| name | admin_level |
| Oregon | 4 |
| Multnomah County | 6 |
| Portland | 8 |
| Portland Downtown | 10 |
這裡需要解釋一下，我們使用的資料集僅限于俄勒岡州，是以它不包含國家/地區等對象。 是以，你不會在結果中看到美國。 其他 OSM 資料導出可能包含更多類型的邊界。
請注意上面 ST_Transform() 函數的使用，以及一些神秘數字，如 4326 和 2269。這就是機關轉換在 PostGIS 中的工作方式。 這些數字稱為空間參考 ID 或 SRID。 SRID 4326 是緯度和經度。 SRID 2269 是一個特定于區域設定的坐标系，稱為“NAD83 / Oregon North (ft)”。 由于地球的形狀有些不規則，而且當你接近兩極時經線會彙聚，是以每個地方都有自己的坐标系。 為了找到距離，我們需要經常将形狀從一個坐标系轉換到另一個坐标系。
4、尋找距離道路最遠的點
現在我們已經掌握了資料的外觀，可以開始工作了。 ST_Distance()、ST_ClosestPoint() 和 ST_Azimuth() 函數适用于任何類型的形狀，而不僅僅是點。 是以，例如，我們可以計算直線上任意一點到最近點的距離。 這将是我們“最遠一條路”計算的核心。 要找到最近的高速公路及其方向，我們可以運作如下查詢：
SELECT osm_id, name, highway, REF,
ST_AsText(ST_Transform(ST_ClosestPoint(way, 'SRID=900913;POINT(-122.679 45.519)'::geometry), 4326)) AS closest_point,
ST_Distance(
way_local,
ST_Transform(ST_GeomFromText('POINT(-122.679 45.519)',4326), 2285)
) AS distance,
degrees(ST_Azimuth(
ST_GeomFromText('POINT(-122.679 45.519)', 4326),
ST_ClosestPoint(ST_Transform(way, 4326), ST_GeomFromText('POINT(-122.679 45.519)', 4326))
)) AS azimuth
FROM roads
WHERE
highway IN ('motorway')
ORDER BY distance
LIMIT 1 現在我們知道與我們的點相關的最近道路在哪裡，我們隻需要對 AI 進行程式設計。
對于這個項目，我們編寫了幾個非常簡單的 AI，它們都遵循類似的計劃：
- 查找到最近道路的距離和方向
- 将點沿相反方向移動一小段距離（例如 1000 英尺）
- 檢查新點是否仍在俄勒岡州内
- 重新計算新點的距離和方向
- 重複這個過程幾次疊代
同樣，它不是太複雜。 如果你将起點放在高速公路環路内，例如環繞波特蘭市中心的 I-5/I-405 環路，它甚至不會聰明到知道它位于封閉形狀内。 是以，該算法也無法“逃脫”封閉形狀以找到其外部的點。 如果離高速公路最遠的點實際上在高速公路的另一邊，則該算法将找不到它。 它所能做的就是找到局部最大距離處的點。 解決這個問題将成為未來AI研究的主題。
5、結束語
一旦我們完成了簡單的 AI，我就有了一個更大的目标。 我想看看我是否能完成我最初的想法，即計算俄勒岡州最偏遠的點。 而且，為了好玩，華盛頓也一樣。
我的計劃是在全州建立一個點網格，間隔大約 2 英裡。 然後我會計算從每個點到任何類型的最近道路的距離。 幾何數運算的計算量很大，是以需要幾個小時才能計算出來。
一旦我确定了最偏遠點的幾個候選點，我就可以使用上述 AI 的修改版本。 我将其程式設計為遠離資料庫中的任何道路，甚至是未命名的四輪驅動軌道。
結果證明這是相當成功的。 我在俄勒岡州東北部（位于美麗的瓦洛厄山脈）确定了距離最近道路 7.5 英裡的一個點，以及華盛頓奧林匹克國家公園距離最近道路 13.6 英裡的一個點。 （我沒有公布确切的位置。我不想冒險讓這些原始地方被遊客淹沒。）
為了驗證我的發現，我在 Project Remote 網站上找到了一張地圖，上面用紅點表示他們在每個州确定的遠端點。 他們的圖像非常縮小，故意不顯示确切位置，但他們在俄勒岡州的東北角和華盛頓州的西北角确實有一個紅點，在我确定的點附近。
我希望你和我一樣覺得這很有趣。
原文連結：http://www.bimant.com/blog/spatial-analysis-with-postgresql-n-openstreetmap/