徐艇偉,張彥
(1.重慶數字城市科技有限公司,重慶 400020; 2.重慶市地理信息云服務企業工程技術研究中心,重慶 400020)
地理國情普查通過對質量元素評分來完成數據的質量檢查。其中,表征質量的折刺檢查可通過電腦程序檢查來完成[1]。地理國情普查的主要數據分為地表覆蓋分類和地理國情要素[2],數據類型包括點狀數據、線狀數據和面狀數據,而面狀數據是最復雜的數據類型,這意味著其折刺檢查工作的難度也最大。[3]通過ArcGIS中折線夾角計算工具完成面狀數據的折刺檢查,存在操作過程復雜、局限性大和效率低下等問題。
隨著現代計算機科學的迅猛發展,計算機在測繪、地理信息系統及遙感信息系統等地理科學行業中的作用越來越大。運用計算機語言,編寫程序來完成檢查是一個解決面狀數據折刺檢查問題的較好途徑。本文通過建立面狀數據的幾何模型,利用計算機的迭代運算和向量夾角計算法來自動獲取每個夾角的角度值,提出一套解決折刺檢查問題的方案。
矢量面狀數據由若干個節點構成,連續3個節點形成一個夾角,當該夾角的角度小于標準閾值時,稱之為折刺。尋找折刺的首要任務是實現面狀數據每個夾角值的自動獲取。
面狀數據可以轉化為點狀數據,通過獲取數據各個節點的位置,利用向量夾角計算法可得到每個節點所對應的夾角度數。向量夾角計算法的具體內容如下:

圖1 三個節點相對位置
假設平面坐標系中的3個節點分別為P0、P1、P2,它們的排列順序如圖1所示,α表示3點連線的夾角,3個點的坐標分別為P0=(x0,y0),P1=(x1,y1),P2=(x2,y2),向量夾角計算公式如下:

通過反三角函數即可獲得向量夾角[4],根據地理國情普查質量評定標準,設定最終夾角小于設定閾值的作為折刺結果導出。
矢量面狀數據折刺的檢查是通過將面狀數據轉化為點狀數據,獲取點的坐標,代入向量夾角計算公式中,得到每個節點夾角的角度值,最后將角度值與設定的標準閾值進行比較,導出小于閾值的角所對應的3個節點的坐標,通過ArcGIS的點生成線工具生成最終的錯誤列表,技術流程圖如圖2所示。

圖2 技術流程圖
以上技術流程的第二步需要將面狀數據轉換為點狀數據,第六步需要生成最終的結果,具體的點狀數據屬性表和結果圖層的屬性表設計如表1和表2所示。

點狀數據屬性結構表 表1

問題點圖層要素屬性結構表 表2
該步驟是將面狀數據中所有節點提取并轉換為點狀數據,通過調用Python的ArcPy模塊中的FeatureVerticesToPoints_management函數可以完成,[5]該功能運算后效果如圖3所示。

圖3 面狀數據轉為點狀數據
(1)獲取每個點的坐標
該步驟通過調用Python的ArcPy模塊中AddXY_management函數,將上一步驟中點狀數據的x和y值計算并記錄到點位屬性列表中。
(2)獲取每個點的前一點和后一點
第一步,獲取每個面起點和終點的信息(起點和終點重合)。由于每個面點狀數據的唯一標識碼從起點到終點是依次遞增的,所以程序通過SQL語言獲取每個面狀數據中包含點狀數據的最小唯一標識碼和最大唯一標識碼。
第二步,根據唯一標識碼獲取每個點的前一點和后一點的信息,每個點前一點的唯一標識碼為該點唯一標識碼減1,每個點后一點的唯一標識碼為該點唯一標識碼加1,由于起點和終點重合的緣故,所以起點和終點的前一點和后一點均是該面的倒數第二個點和整數第二個點。
該步驟為本檢查方法的核心步驟,具體算法流程如圖4所示。

圖4 算法流程圖
(1)利用Python的數組List功能和Arcpy模塊的Cursor功能,獲取面的個數得到N值,關鍵代碼如下:
polygons[]
with arcpy.da.SearchCursor(fc,["OID@","ORIG_FID"]) as cursor:#檢索每個面狀唯一標識碼
for row in cursor:
polygons.append(row[1])
del cursor
N=len(polygons) #獲取面的個數
(2)獲取每個面起點和終點的唯一標識碼,關鍵代碼如下:
for m in range(1,N):
with arcpy.da.SearchCursor(fc,["ORIG_FID","MIN_OBJECTID","MAX_OBJECTID"],"ORIG_FID="+str(m)) as cursor: #檢索每個面狀數據的唯一標識碼、所含點狀數據的最小標識碼和最大標識碼
for row in cursor:
sta=row[1] #獲取起點ID
end=row[2] #獲取終點ID
……
(3)獲取當前面起點、第2點、第3點……第n點……終點的x值、y值、前一點唯一標識碼、后一點唯一標識碼,關鍵代碼如下:
with arcpy.da.SearchCursor(fc,["OID@","ORIG_FID","POINT_X","POINT_Y","Prev","Next"],"OBJECTID="+str(sta)) as cursor: #獲取起點的信息
for row in cursor:
x1=row[2] y1=row[3] p1=row[4] n1=row[5]
with arcpy.da.SearchCursor(fc,["OID@","ORIG_FID","POINT_X","POINT_Y","Prev","Next"],"OBJECTID="+str(n1)) as cursor: #獲取第2點的信息
for row in cursor:
cid=row[0] x2=row[2] y2=row[3] cp=row[4] cn=row[5]
while True:
with arcpy.da.SearchCursor(fc,["OID@","POINT_X","POINT_Y","Next"],"OBJECTID="+str(cn)) as cursor:#獲取第n點的信息(n≠1,2)
for row in cursor:
nid=row[0] nx=row[1] ny=row[2] cn=row[3][6]
(4)求當前面中第n個夾角值的余弦值,并通過反三角函數求得第n個夾角的角度值。構成該夾角的三個點分別為第n-1點,第n點,第n+1點,坐標值分別為Pn-1=((n-1)x,(n-1)y),Pn=(nx,ny),Pn+1=((n+1)x,(n+1)y),關鍵代碼如下:
(((n-1)x-nx)*((n+1)x-nx)+((n-1)y-ny)*((n+1)y-ny))/(((n-1)x-nx)**2+((n-1)y-ny)**2)**0.5/(((n+1)x-nx)**2+((n+1)y-ny)**2)**0.5
(5)通過條件語句對第n個夾角余弦值與標準值的余弦值進行比較,將構成超限夾角的三個點狀數據導出到點問題圖層中,形成初步檢查結果。
為了方便檢查結果顯示,將點問題圖層轉化為線狀數據,通過Python中ArcPy模塊的點生成線函數可以完成。
本文采用第一次地理國情普查重慶市永川區的地表覆蓋分類數據來驗證面狀數據折刺檢查方法。永川區地表覆蓋分類數據為面狀數據,試驗區域包含面狀數據共3 857個,包含節點 125 273個。
(1)將普查數據中的地表覆蓋分類數據轉換為點狀數據,如圖5所示。

圖5面狀數據轉化為點狀數據
(2)獲取點位坐標,計算3點連線的夾角,生成結果如表3所示。

向量夾角計算結果表 表3
表中記錄第1824、1825、1826三個點和第7896、7897、7898三個點所構成的兩個夾角值分別約為13.68°和6.07°,根據地理國情普查地表覆蓋分類檢查標準中面狀數據折刺角度不得小于15°的規定,這六個點涉及的兩個夾角均為異常折刺,作為問題點位導出。
(3)導出超限夾角,生成最終結果圖層,如圖6所示。

圖6 問題線圖層
上述方法與ArcGIS中的折線夾角計算方法進行對比,有三大優點:
(1)適用性更廣,ArcGIS的折點夾角計算只適用于線型數據,而上述方法既可以計算線型數據也可以計算面狀數據;
(2)易用性更高,ArcGIS中的折點夾角計算方法需要將數據轉換為點狀數據后,建立COGO字段,通過COGO工具獲取線段的方向角,從而得到折線夾角[7]。而利用本文論述的方法,只需要導入待檢測的數據,即可運算出檢查結果,不需要介入更多的人工操作。
(3)效率更高,如果使用ArcGIS的折點夾角計算方法,需要耗費的時間遠比使用上述方法要長。
矢量面狀數據折刺自動檢查方法的優勢是從根本上解決了基礎地理數據質量檢查的角度值異常問題,也為其他基礎地理數據幾何問題的排查提供了解決思路。該方法不僅可以運用在面狀數據的折刺檢查中,同樣可運用在線狀數據的折刺檢查中。但同時也存在弊端,隨著運算量增大,運算速度降低更快,當數據較大時,可通過劃分查詢區域、精簡查詢步驟以及利用已有查詢結果等方法提升運算效率。計算機語言中,除了PYTHON語言,也可通過其他編程語言和平臺來進行自動檢查方面的開發和研究,借助計算機來完成質量檢查工作將成為未來地理信息行業質檢領域的一大趨勢。