有效點在多邊形連接通過PySpark和BNG地理空間索引

(不)隱藏的成本?

在討論加入關係複雜性,我們取得了一個監督。傳統的複雜性假設每一對決議的固定成本,也就是說,到達一個結論的成本每一對匹配或不匹配的記錄在連接操作,我們將調用O(加入)。加入的真實成本O (n *米)* O(加入)。在傳統的等價類的關係,我們隻是在考慮是否加入鍵左邊匹配連接鍵在右邊,我們假設O(加入)是O (1)或者簡單地說,比較的成本是一個算術運算,它是恒定的。這並非總是如此;例如,加入一個字符串比較是更昂貴的比兩個整數之間的等價性。

但是皮普,相對的成本如何?使用最廣泛的算法回答是皮普文中算法。該算法的複雜性O (v),v是多邊形的頂點的數量問題。該算法適用於凸和非凸形狀,和保持相同的複雜性在這兩種情況下。

增加的比較成本成本模型使我們總立方形式的複雜性。如果我們更換O(加入)與O (v)v是頂點的平均數量,我們總的複雜性O (n *米)* O (v)。這既昂貴又耗時的!

更聰明地工作,而不是更刻苦!

我們可以做得更好O (n *米)* O (v)。我們可以使用火花來幫助我們擊敗了笛卡兒積的複雜性。火花利用散列連接。根據連接謂詞,火花可以執行下列之一連接策略:

• 廣播散列連接的複雜性O (max (n, m)) * O(加入)
• 隨機散列連接(類似於優雅散列連接)的複雜性O (n + m) * O(加入)
• 洗牌分類合並加入的複雜性O (n * log (n) + m *日誌(m)) * O(加入)
• 廣播嵌套循環連接(笛卡爾連接)和複雜性O (n *米)* O(加入)

神奇的!我們可以使用火花,我們將避免最昂貴的結果,我們不能?不!不幸的是,火花將默認為PIP笛卡爾連接連接。為什麼?在皮普有別於傳統equi-based連接是基於一般的關係。這些連接通常稱為θ加入類。這些通常是難以執行,需要最終用戶幫助係統。當我們從一個不利的位置,我們仍然可以實現所需的性能。

空間索引pseudo-equivalence (PIP)

有辦法讓皮普是一個等價關係?嚴格地說不,然而,在實踐中,我們可以讓皮普方法效率的一個等價關係,如果我們采用空間索引技術。

空間指標有效地幫助我們指數坐標空間的邏輯分組中接近另一個幾何圖形表示空間。我們實現這一獨特的關聯點的坐標係統索引ID。這些係統允許我們代表參考空間不同層次的細節,或簡單,不同的分辨率。此外,地理空間索引係統分層係統;這意味著有一個定義良好的親子關係指數在不同級別的表示。

這將如何幫助我們嗎?如果我們分配給每個幾何圖形它所屬的指數,我們可以使用索引ID來索引ID等效為一個等價關係代理。我們將執行脈衝(或任何其他幾何投影關係)隻屬於同一指標的幾何圖形。

重要的是要注意,雖然點幾何圖形屬於一個且隻有一個指數,所有其他的幾何類型,包括線和多邊形,可能跨度超過一組的指數。這意味著成本的解決脈衝通過索引空間關係O (k) * O (v),k指數用來表示數量的幾何和v是頂點的幾何的數量。這表明我們每個的價格比較增加爆炸記錄複雜的幾何圖形為多個索引記錄攜帶相同的幾何。

為什麼這是一個明智的選擇嗎?當我們正在增加的價格比較單一的幾何圖形,我們避免了一個完整的笛卡兒積,還是在大型地理空間連接。稍後我們將詳細顯示,索引ID索引ID連接將使我們能夠跳過大量不必要的比較。

最後,數據源包含複雜幾何形狀不發展盡快做逐點數據源。複雜的幾何圖形通常代表地區,感興趣的領域,建築,等等,這些概念有相當穩定的時間表,對象隨時間變化的變化很少,和對象變化相對較少。這意味著,當我們做的花費額外的時間進行預處理複雜的幾何圖形,對大部分人來說,這預處理是一個一次性的事件。這種方法仍然適用甚至頻繁更新數據;我們可以跳過,當加入的數據量通過索引ID來索引ID大於增加的行數的關係用來表示一個幾何。

BNG指標體係

BNG是地方坐標參考係統(CRS) (EPSG: 27700)成立於1936年,為國家映射,包括英國。與全球CRS, BNG已經安裝和塑造英國、大陸的坐標投影到平麵,普通方格網的原點(0,0)的西南錫利群島的。

在網格範圍內,地理網格引用(或指標)是用來確定網格方塊在不同的決議表示在米可以翻譯和BNG以東(x)和北航(y)坐標。鑒於網格原點的位置,以東和以北值值總是正的。BNG作為主要參考係統下的所有操作係統位置數據捕獲他們的國家公共任務映射,因此,被廣泛采用的公共和私人用戶操作係統數據操作在英國。

每一個網格都可以表示為一個多邊形幾何體,每條邊的長度等於參考網格的分辨率。這使得BNG更容易起點地理空間數據分區策略。我們的起點是一個正方形作為構建塊,和它會讓很多開始考慮簡單的同時不丟失的推廣方法。

按照慣例,BNG網格引用表示為字符串,使用字母和西南角的一個給定的網格坐標引用一個特定的解決方案。任何參考的前兩個字符是字母(前綴)(例如,TQ)識別91網格的廣場之一測量100.000米(100公裏)。隻有55 91 100公裏網格方塊介紹一些在英國大陸。餘下的這些方塊落入英國水域。

引用識別更細粒度的網格分辨率低於100公裏將會有額外的x和y整數值後附加兩個字母定位子網格是在父網格層次結構。兒童廣場從左下編號從0到9(西南)的角落,在一個東風(x)和北方(y)方向。

為什麼BNG ?

同時有替代全局索引係統,我們可以采用這項工作,我們選擇使用BNG因為:

• BNG係統是本機操作係統的地理空間數據收集,與幾乎所有的操作係統數據引用對BNG CRS (EPSG: 27700)。這包括操作係統空中影像瓦片等柵格數據集,如數字地形模型(dtm)和數字表麵模型(dsm)。
• BNG使高效檢索的使用和主機托管的矢量和柵格數據分析,包括剪切或屏蔽的柵格數據推導培訓深度學習應用程序補丁,作為一個例子。
• 使用BNG避免了昂貴的世界大地轉換係統(wgs - 84) (1984EPSG: 4326)或1989年歐洲地麵參考係統(ETRS89) (EPSG: 4258)crs通過OSTN15網格轉換。不同crs實現他們的地球模型使用不同的參數,和全球係統(例如,WGS84)將顯示一個偏移量相比,一個本地係統(例如,BNG)。這種轉換的真實成本反映在操作係統這一事實OSTN15出版,一個包含大約15 mb修正文件。175萬參數采用衛星坐標和BNG坐標之間的變換準確。

由於GB-local操作係統試圖解決的問題的本質,BNG是一種自然選擇。在更多的全局上下文的情況下,我們應該轉換我們的關注H3或S2成為全球選擇更合適。

BNG空間劃分策略

空間劃分策略定義了一個方法將地理空間數據分割到非重疊區域。BNG網格方塊在不同的決議在英國提供非重疊區域。通過檢索BNG指標,涵蓋幾何圖形我們可以使用索引屬性搭配加入關鍵行然後隻測試一個空間謂詞在這些集中的行(例如,幾何相交的幾何B還是幾何包含幾何B)。

這是非常重要的!把原始數據分解為地理空間數據集中的部分“高度平行”,和使我們的問題,因此,非常適合火花/ PySpark。我們可以將不同的數據塊發送給不同的機器,隻比較當地的部分數據,可能會加入一個到另一個地方。有小點檢查建築在倫敦包含在曼徹斯特的一個地址。地理空間指標是我們的方式來傳達這種直覺的機器。

基線

我們使用Python和PySpark給生活帶來我們的解決方案。操作係統提供的邏輯轉換提供的坐標經常和北航獨特BNG索引ID。最後,為了確保一個公正的輸出,我們使用一個隨機點和多邊形的隨機數據集的數據集;1000萬點散落在GB,境內100萬多邊形被分散在同樣的方式。生成一組多邊形數據,我們已經加載一套GeoJSON dataframe火花,我們使用一個隨機函數與一個生成器函數(爆炸)來生成一個公正的數據集。由於隨機性引入數據,應當期望點和多邊形之間的關係是多對多的。

基線算法用於我們的考慮是天真的加入,會導致非θ加入。這種方法將執行時間,被評價為一個嵌套循環聯接播放。

廣播嵌套循環聯接運行非常緩慢。這是事實,原因是評估同樣笛卡爾連接。每個point-polygon對評估對皮普關係加入之前解決。結果是,我們需要加入十億比較10萬點到1萬多邊形。注意,這兩個數據集是大到足以被稱為大數據。

我們使用MLflow進行一係列的天真加入評估基線性能我們試圖超越。天真的方法,最大的加入我們能夠成功執行1萬點到10萬多邊形。數據量的進一步增加導致我們的火花工作沒有產生所需的輸出。這些失敗是由於實現工作負載的性質我們試圖運行。

讓我們框架的問題

如果我們代表我們所有的幾何圖形,無論其形狀,與相應的BNG-aligned邊界框嗎?邊界框是一個矩形多邊形可以符合原文的整體幾何。如果我們說邊界框表示為一組BNG指數在給定的決議,覆蓋同一區域。

現在我們可以執行連接通過一個更優化的θ加入。我們隻檢查一個點是否在多邊形通過皮普關係如果一個點落入一個BNG指數用於表示多邊形。這樣可以減少我們的加入工作由多個數量級。

為了生產說BNG指數,我們使用下麵的代碼;請注意,bng_to_geom, coords_to_bng bng_get_resolution功能不提供這個博客。

從shapely.geometry進口盒子#輔助函數來檢索第一個鄰國# BNG指數的細胞defnext_horizontal(bng_index,決議):x, y = bng_to_geom (bng_index)返回coords_to_bng (x +決議,y,分辨率)#輔助函數來檢索第一個鄰國# BNG指數細胞的底部defnext_vertical(bng_index,決議):x, y = bng_to_geom (bng_index)返回coords_to_bng (x, y-resolution分辨率)#填充函數,表示輸入幾何的一組指標#對應區域的邊界框的幾何表示defbng_polyfil(多邊形,決議):(x1, y1, x2, y2) = polygon.boundsbounding_box =盒(* polygon.bounds)lower_left = coords_to_bng (x1, y2,分辨率)隊列= [lower_left]結果=集()參觀了=集()而隊列:指數= queue.pop ()index_geom = shapely.wkt.loads (bng_to_geom_grid(指數,“WKT”))十字路口= bounding_box.intersects (index_geom)如果十字路口:result.add(索引)n_h = next_horizontal(指數、分辨率)如果n_h不在訪問:queue.append (n_h)n_v = next_vertical(指數、分辨率)如果n_v不在訪問:queue.append (n_v)visited.add(索引)訪問= []返回結果

這段代碼確保我們可以無損地代表任何形狀。我們使用BNG指數候選人相交關係,避免盲點的原始幾何表示。注意,更有效的實現是可能通過使用包含關係和質心點;這種方法是唯一可行的是否可以接受的假陽性和假陰性。我們假設的存在bng_to_geom函數,給定一個BNG索引ID可以產生一個幾何表示,bng_get_resolution函數,給定一個BNG決定選擇的分辨率和索引IDcoords_to_bng返回一個函數,考慮到坐標BNG索引ID。

我們運行我們的多邊形邊界框表示BNG指標體係的不同分辨率和不同大小的數據集。注意,運行這個過程失敗一直低於100號決議。在這些輸出決議米表示。一致的失敗的原因在代表的決議可以找到低於100米;一些多邊形(由於隨機性質)比其他的大得多,而一些多邊形是由一組12個指標,其他多邊形可以表示為成千上萬的指標,這可能會導致一個大的差距在計算和內存需求在分區之間的火花工作產生這些數據。

我們省略了因為這個點的基準數據集轉換是一個相對簡單的操作,不會帶來任何新行;隻添加一個列,不同的分辨率不影響執行時間。

雙方的加入是用相應的BNG表示表示,所有我們要做的就是執行調整後加入邏輯:

@udf (“布爾”)defpip_filter(poly_wkt、point_x point_y):從有條理的進口wkt從有條理的進口幾何多邊形= wkt.loads (poly_wkt)點=幾何。點(point_x point_y)返回polygon.contains(點)defrun_bounding_box_join(polygons_path, points_path):多邊形= spark.read。格式(“δ”).load (polygons_path)多邊形= polygons.select (F.col (“id”),F.col (“wkt_polygon”),F.explode (F.col (“bng_set”).alias (“bng”))點= spark.read。格式(“δ”).load (points_path)返回polygons.join (點,= (“bng”),如何=“內心”)。(pip_filter (“wkt_polygon”,“經常”,“北航”))#連接數據集上運行一個動作來評估加入運行時run_bounding_box_join (polygons_path points_path) .count ()

這些修改我們的代碼導致了不同的火花執行計劃。火花現在能夠首次運行基於BNG指數排序合並連接ID和大大減少的總數比較。此外,每一對比較string-to-string比較遠短於一個脈衝的關係。第一階段將生成所有加入組候選人。我們將執行PIP關係測試組候選人解決最終的輸出。這種方法可以確保我們限製的次數PIP操作運行。

執行計劃,我們可以看到,火花是執行不同的操作相比,幼稚的方法。最值得注意的是,火花現在執行排序合並連接,而不是廣播嵌套循環連接,這是讓很多效率。我們現在執行大約1.86億PIP操作而不是十億年。這就允許我們運行更大的連接有更好的響應時間,同時避免任何突發故障,我們經曆了天真的方法。

這個簡單但有效的優化使我們運行一個PIP加入1000萬點至100萬多邊形在大約2500秒。如果我們比較基線執行時間,最大的加入我們能夠成功執行1萬點到10萬多邊形,甚至加入需要大約1500秒在相同的硬件上。

分而治之

運行數據集在幾百萬行域之間的連接是偉大的;然而,我們最大的基準加入了將近45分鍾(2500秒)。在世界上我們想要運行特設在大量地理空間數據分析,這些執行時間太緩慢。

我們需要進一步優化的方法。第一個候選人優化我們的邊界框表示。如果我們通過邊界框代表多邊形,我們有太多的假陽性指標,即。指數,與原來的幾何不重疊。

優化部分的代碼的方法是簡單地使用在我們polyfill相交函數調用方法的原始幾何。

defk_ring(bng_index):x, y = bng_to_geom (bng_index)增量= bng_get_resolution (bng_index)鄰居= [[x-increment y +增量],[x, y +增量],[x +增量,y +增量)[x-increment y], [x +增量,y],[x-increment, y-increment], [x, y-increment], [x +增量,y-increment]]鄰居= [coords_to_bng(我0),我1),增量)為我在鄰居)返回鄰居defbng_polyfil(多邊形,決議):從shapely.geometry進口盒子開始= get_starting_point(多邊形,分辨率)隊列= k_ring(開始)結果=集()參觀了=集()而隊列:指數= queue.pop ()如果polygon.intersects (shapely.wkt.loads (bng_to_geom_grid(指數,“WKT”))):result.add(索引)為n在k_ring(指數):如果n不在參觀了和n不在隊列:queue.append (n)visited.add(索引)訪問= []返回結果

這種優化,利用成本增加相交調用,將導致更小的結果指標集和將使我們連接運行得更快,由於小加入表麵

第二優化我們可以使用表示分解為兩套指標。並不是所有的指標都是平等的在我們的表示。指數的接觸邊界多邊形需要脈衝過濾索引索引後加入。指數,請勿觸摸多邊形的邊界和屬於表示不需要任何額外的過濾。任何時候,落入這樣一個指數絕對屬於多邊形,在這種情況下,我們可以跳過PIP操作。

我們可以實現的第三個和最後一個優化是馬賽克的方法。而不是將完整的原始幾何與屬於的指標集的每個索引聯係多邊形邊界(邊境),我們隻能跟蹤感興趣的部分。如果我們相交的幾何表示和多邊形指數問題,我們得到的地方表示多邊形;隻有這部分的原始多邊形的麵積指數的相關問題。我們稱這些作品為多邊形芯片。

從優化角度多邊形芯片有兩個目的。首先,他們大大提高PIP過濾器的效率發生index-to-index加入後執行。這是由於這樣的事實,射線跟蹤算法運行在O (v)的複雜性和個體芯片平均少一個數量級比原始幾何頂點。其次,芯片的代表性遠小於原來的幾何形狀,因此,我們拖著更少的數據作為洗牌階段的一部分在我們的排序合並連接階段。

把所有這些一起收益率以下代碼:

defadd_children(隊列、訪問索引):為n在k_ring(指數):如果n不在參觀了和n不在隊列:queue.append (n)返回隊列defbng_polyfil(多邊形,決議):開始= get_starting_point(多邊形,分辨率)隊列= k_ring(開始)結果=集()參觀了=集()而隊列:指數= queue.pop ()index_geom = shapely.wkt.loads (bng_to_geom_grid(指數,“WKT”))十字路口= polygon.intersection (index_geom)如果intersection.equals (index_geom):result.add(指數,假,“多邊形空”))隊列= add_children(隊列、訪問、索引)elif“空”不在intersection.to_wkt ():result.add(指數,真正的intersection.to_wkt ()))隊列= add_children(隊列、訪問、索引)visited.add(索引)訪問= []返回結果

這段代碼非常類似於原始的邊界框的方法,我們隻做了一些微小的變化,確保我們不是複製部分代碼;因此,我們孤立add_children輔助方法。

執行相同的數據生成基準我們為我們所做的邊界框多邊形表示。我們發現一個共同點與原方法分辨率低於100米導致多邊形的代表。然而,在這種情況下,我們能夠生成數據多達10萬個多邊形10米的分辨率,獲得運行時的數據生成過程太緩慢,被認為是生產工作負載。

在100的分辨率,我們有一些非常不錯的效果;花了大約600秒的數據集生成並寫出100萬多邊形。作為參考,花了大約300秒做同樣的邊界框的方法。邊界框是一個簡單的過程,我們添加一些處理時間在數據準備階段。我們可以證明這些投資嗎?

馬賽克是相當(快!)

我們已經運行相同的基準脈衝連接使用我們的馬賽克數據。我們改編加入邏輯略以確保我們的邊境設置和核心指標都是利用正確和最有效的方式。

def run_polyfill_chipping_join (polygons_path points_path):多邊形=spark.read.format(“δ”).load (polygons_path)多邊形=polygons.select (F.col (" id "),F.explode (F.col (“bng_set”)) .alias (“bng”))。選擇(F.col (" id "),F.col (“bng *”。))正確的=spark.read.format(“δ”).load (right_path)返回polygons.join (正確的,在=[" bng "),如何=“內心”)。在哪裏(~F.col (“is_dirty”)|pip_filter (“wkt_chip”、“經常”、“北航”))#運行一個行動來執行的加入run_polyfill_chipping_join (polygons_path points_path)。數()

is_dirty列是由我們polyfill方法引入。任何觸及邊境的原始幾何指數將標記為髒(例如,is_dirty = True)。這些指標要求後過濾,以便正確地確定任意點,落入指數比較中包含的幾何表示。這是至關重要的is_dirty過濾發生之前pip_fiter邏輯運算符調用,因為在火花短路能力;如果第一部分的邏輯表達式是正確的,第二部分將不會執行。

這段代碼將產生一個更有效的執行計劃的火花。由於索引空間更好的表現,我們加入表麵要小得多。此外,我們post-filters受益於2集表示和馬賽克幾何圖形的分割。

我們終於可以量化我們的努力。PIP類型加入1000萬點和100萬個多邊形之間通過我們新的馬賽克方法已經在37秒內執行。把這個到上下文中,等效邊界框連接在同一索引決議執行2549秒。這導致69 x改善運行時。

這種改進單純側重於運行時服務。如果我們包括準備時間,600秒的鑲嵌方法,317秒的邊界框的方法,我們有總調整性能提高4.5倍。

這些改進的總潛力很大程度上取決於你有多經常更新你的幾何數據和多久你查詢它。

通用的方法

在這篇文章中,我們主要集中在點在多邊形(PIP)連接使用英國國家電網(BNG)作為參考指標體係。但是,比這更一般的方法。相同的優化可以適應任何層次地理空間係統。所不同的是,芯片的形狀和可用的決議。此外,相同的優化可以幫助你擴大θ之間連接兩個複雜的幾何圖形,如大量多邊形交叉連接。

我們依然關注PySpark第一種方法,有意識地避免引入任何第三方框架。我們相信,消耗我們的解決方案,確保較低的障礙主要是Python用戶量身定製。

幾個創意優化我們的解決方案已經證明,可以達到70倍的性能改進與最小邊界框的方法預處理增加投資。

我們帶來了大規模的PIP加入到執行時間域的秒,我們打開特別針對這些數據的分析功能。