徐艇偉，張彥

(1.重慶數(shù)字城市科技有限公司，重慶 400020； 2.重慶市地理信息云服務(wù)企業(yè)工程技術(shù)研究中心，重慶 400020)

1 引 言

地理國情普查通過對質(zhì)量元素評分來完成數(shù)據(jù)的質(zhì)量檢查。其中，表征質(zhì)量的折刺檢查可通過電腦程序檢查來完成[1]。地理國情普查的主要數(shù)據(jù)分為地表覆蓋分類和地理國情要素[2]，數(shù)據(jù)類型包括點狀數(shù)據(jù)、線狀數(shù)據(jù)和面狀數(shù)據(jù)，而面狀數(shù)據(jù)是最復(fù)雜的數(shù)據(jù)類型，這意味著其折刺檢查工作的難度也最大。[3]通過ArcGIS中折線夾角計算工具完成面狀數(shù)據(jù)的折刺檢查，存在操作過程復(fù)雜、局限性大和效率低下等問題。

隨著現(xiàn)代計算機科學(xué)的迅猛發(fā)展，計算機在測繪、地理信息系統(tǒng)及遙感信息系統(tǒng)等地理科學(xué)行業(yè)中的作用越來越大。運用計算機語言，編寫程序來完成檢查是一個解決面狀數(shù)據(jù)折刺檢查問題的較好途徑。本文通過建立面狀數(shù)據(jù)的幾何模型，利用計算機的迭代運算和向量夾角計算法來自動獲取每個夾角的角度值，提出一套解決折刺檢查問題的方案。

2 方法原理與設(shè)計

2.1 原理分析

矢量面狀數(shù)據(jù)由若干個節(jié)點構(gòu)成，連續(xù)3個節(jié)點形成一個夾角，當(dāng)該夾角的角度小于標(biāo)準(zhǔn)閾值時，稱之為折刺。尋找折刺的首要任務(wù)是實現(xiàn)面狀數(shù)據(jù)每個夾角值的自動獲取。

面狀數(shù)據(jù)可以轉(zhuǎn)化為點狀數(shù)據(jù)，通過獲取數(shù)據(jù)各個節(jié)點的位置，利用向量夾角計算法可得到每個節(jié)點所對應(yīng)的夾角度數(shù)。向量夾角計算法的具體內(nèi)容如下：

圖1 三個節(jié)點相對位置

假設(shè)平面坐標(biāo)系中的3個節(jié)點分別為P0、P1、P2，它們的排列順序如圖1所示，α表示3點連線的夾角，3個點的坐標(biāo)分別為P0=(x0，y0)，P1=(x1，y1)，P2=(x2，y2)，向量夾角計算公式如下：

通過反三角函數(shù)即可獲得向量夾角[4]，根據(jù)地理國情普查質(zhì)量評定標(biāo)準(zhǔn)，設(shè)定最終夾角小于設(shè)定閾值的作為折刺結(jié)果導(dǎo)出。

2.2 技術(shù)設(shè)計

矢量面狀數(shù)據(jù)折刺的檢查是通過將面狀數(shù)據(jù)轉(zhuǎn)化為點狀數(shù)據(jù)，獲取點的坐標(biāo)，代入向量夾角計算公式中，得到每個節(jié)點夾角的角度值，最后將角度值與設(shè)定的標(biāo)準(zhǔn)閾值進行比較，導(dǎo)出小于閾值的角所對應(yīng)的3個節(jié)點的坐標(biāo)，通過ArcGIS的點生成線工具生成最終的錯誤列表，技術(shù)流程圖如圖2所示。

圖2 技術(shù)流程圖

2.3 屬性設(shè)計

以上技術(shù)流程的第二步需要將面狀數(shù)據(jù)轉(zhuǎn)換為點狀數(shù)據(jù)，第六步需要生成最終的結(jié)果，具體的點狀數(shù)據(jù)屬性表和結(jié)果圖層的屬性表設(shè)計如表1和表2所示。

點狀數(shù)據(jù)屬性結(jié)構(gòu)表 表1

問題點圖層要素屬性結(jié)構(gòu)表 表2

3 程序設(shè)計

3.1 面狀數(shù)據(jù)轉(zhuǎn)化為點狀數(shù)據(jù)

該步驟是將面狀數(shù)據(jù)中所有節(jié)點提取并轉(zhuǎn)換為點狀數(shù)據(jù)，通過調(diào)用Python的ArcPy模塊中的FeatureVerticesToPoints_management函數(shù)可以完成，[5]該功能運算后效果如圖3所示。

圖3 面狀數(shù)據(jù)轉(zhuǎn)為點狀數(shù)據(jù)

3.2 獲取點的坐標(biāo)及其前一點和后一點

(1)獲取每個點的坐標(biāo)

該步驟通過調(diào)用Python的ArcPy模塊中AddXY_management函數(shù)，將上一步驟中點狀數(shù)據(jù)的x和y值計算并記錄到點位屬性列表中。

(2)獲取每個點的前一點和后一點

第一步，獲取每個面起點和終點的信息(起點和終點重合)。由于每個面點狀數(shù)據(jù)的唯一標(biāo)識碼從起點到終點是依次遞增的，所以程序通過SQL語言獲取每個面狀數(shù)據(jù)中包含點狀數(shù)據(jù)的最小唯一標(biāo)識碼和最大唯一標(biāo)識碼。

第二步，根據(jù)唯一標(biāo)識碼獲取每個點的前一點和后一點的信息，每個點前一點的唯一標(biāo)識碼為該點唯一標(biāo)識碼減1，每個點后一點的唯一標(biāo)識碼為該點唯一標(biāo)識碼加1，由于起點和終點重合的緣故，所以起點和終點的前一點和后一點均是該面的倒數(shù)第二個點和整數(shù)第二個點。

3.3 計算向量夾角并比較

該步驟為本檢查方法的核心步驟，具體算法流程如圖4所示。

圖4 算法流程圖

(1)利用Python的數(shù)組List功能和Arcpy模塊的Cursor功能，獲取面的個數(shù)得到N值，關(guān)鍵代碼如下：

polygons[]

with arcpy.da.SearchCursor(fc，["OID@"，"ORIG_FID"]) as cursor：#檢索每個面狀唯一標(biāo)識碼

for row in cursor：

polygons.append(row[1])

del cursor

N=len(polygons) #獲取面的個數(shù)

(2)獲取每個面起點和終點的唯一標(biāo)識碼，關(guān)鍵代碼如下：

for m in range(1，N)：

with arcpy.da.SearchCursor(fc，["ORIG_FID"，"MIN_OBJECTID"，"MAX_OBJECTID"]，"ORIG_FID="+str(m)) as cursor： #檢索每個面狀數(shù)據(jù)的唯一標(biāo)識碼、所含點狀數(shù)據(jù)的最小標(biāo)識碼和最大標(biāo)識碼

for row in cursor：

sta=row[1] #獲取起點ID

end=row[2] #獲取終點ID

……

(3)獲取當(dāng)前面起點、第2點、第3點……第n點……終點的x值、y值、前一點唯一標(biāo)識碼、后一點唯一標(biāo)識碼，關(guān)鍵代碼如下：

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)求當(dāng)前面中第n個夾角值的余弦值，并通過反三角函數(shù)求得第n個夾角的角度值。構(gòu)成該夾角的三個點分別為第n-1點，第n點，第n+1點，坐標(biāo)值分別為Pn-1=((n-1)x，(n-1)y)，Pn=(nx，ny)，Pn+1=((n+1)x，(n+1)y)，關(guān)鍵代碼如下：

(((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個夾角余弦值與標(biāo)準(zhǔn)值的余弦值進行比較，將構(gòu)成超限夾角的三個點狀數(shù)據(jù)導(dǎo)出到點問題圖層中，形成初步檢查結(jié)果。

3.4 生成問題線圖層

為了方便檢查結(jié)果顯示，將點問題圖層轉(zhuǎn)化為線狀數(shù)據(jù)，通過Python中ArcPy模塊的點生成線函數(shù)可以完成。

4 試驗結(jié)果與分析

4.1 試驗數(shù)據(jù)

本文采用第一次地理國情普查重慶市永川區(qū)的地表覆蓋分類數(shù)據(jù)來驗證面狀數(shù)據(jù)折刺檢查方法。永川區(qū)地表覆蓋分類數(shù)據(jù)為面狀數(shù)據(jù)，試驗區(qū)域包含面狀數(shù)據(jù)共3 857個，包含節(jié)點 125 273個。

4.2 試驗過程

(1)將普查數(shù)據(jù)中的地表覆蓋分類數(shù)據(jù)轉(zhuǎn)換為點狀數(shù)據(jù)，如圖5所示。

圖5面狀數(shù)據(jù)轉(zhuǎn)化為點狀數(shù)據(jù)

(2)獲取點位坐標(biāo)，計算3點連線的夾角，生成結(jié)果如表3所示。

向量夾角計算結(jié)果表 表3

表中記錄第1824、1825、1826三個點和第7896、7897、7898三個點所構(gòu)成的兩個夾角值分別約為13.68°和6.07°，根據(jù)地理國情普查地表覆蓋分類檢查標(biāo)準(zhǔn)中面狀數(shù)據(jù)折刺角度不得小于15°的規(guī)定，這六個點涉及的兩個夾角均為異常折刺，作為問題點位導(dǎo)出。

(3)導(dǎo)出超限夾角，生成最終結(jié)果圖層，如圖6所示。

圖6 問題線圖層

4.3 試驗分析

上述方法與ArcGIS中的折線夾角計算方法進行對比，有三大優(yōu)點：

(1)適用性更廣，ArcGIS的折點夾角計算只適用于線型數(shù)據(jù)，而上述方法既可以計算線型數(shù)據(jù)也可以計算面狀數(shù)據(jù)；

(2)易用性更高，ArcGIS中的折點夾角計算方法需要將數(shù)據(jù)轉(zhuǎn)換為點狀數(shù)據(jù)后，建立COGO字段，通過COGO工具獲取線段的方向角，從而得到折線夾角[7]。而利用本文論述的方法，只需要導(dǎo)入待檢測的數(shù)據(jù)，即可運算出檢查結(jié)果，不需要介入更多的人工操作。

(3)效率更高，如果使用ArcGIS的折點夾角計算方法，需要耗費的時間遠(yuǎn)比使用上述方法要長。

5 總 結(jié)

矢量面狀數(shù)據(jù)折刺自動檢查方法的優(yōu)勢是從根本上解決了基礎(chǔ)地理數(shù)據(jù)質(zhì)量檢查的角度值異常問題，也為其他基礎(chǔ)地理數(shù)據(jù)幾何問題的排查提供了解決思路。該方法不僅可以運用在面狀數(shù)據(jù)的折刺檢查中，同樣可運用在線狀數(shù)據(jù)的折刺檢查中。但同時也存在弊端，隨著運算量增大，運算速度降低更快，當(dāng)數(shù)據(jù)較大時，可通過劃分查詢區(qū)域、精簡查詢步驟以及利用已有查詢結(jié)果等方法提升運算效率。計算機語言中，除了PYTHON語言，也可通過其他編程語言和平臺來進行自動檢查方面的開發(fā)和研究，借助計算機來完成質(zhì)量檢查工作將成為未來地理信息行業(yè)質(zhì)檢領(lǐng)域的一大趨勢。