徐艇偉，張彥

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

1 引 言

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

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

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

2.1 原理分析

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

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

圖1 三個(gè)節(jié)點(diǎn)相對(duì)位置

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

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

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

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

圖2 技術(shù)流程圖

2.3 屬性設(shè)計(jì)

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

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

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

3 程序設(shè)計(jì)

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

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

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

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

(1)獲取每個(gè)點(diǎn)的坐標(biāo)

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

(2)獲取每個(gè)點(diǎn)的前一點(diǎn)和后一點(diǎn)

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

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

3.3 計(jì)算向量夾角并比較

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

圖4 算法流程圖

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

polygons[]

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

for row in cursor：

polygons.append(row[1])

del cursor

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

(2)獲取每個(gè)面起點(diǎn)和終點(diǎn)的唯一標(biāo)識(shí)碼，關(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： #檢索每個(gè)面狀數(shù)據(jù)的唯一標(biāo)識(shí)碼、所含點(diǎn)狀數(shù)據(jù)的最小標(biāo)識(shí)碼和最大標(biāo)識(shí)碼

for row in cursor：

sta=row[1] #獲取起點(diǎn)ID

end=row[2] #獲取終點(diǎn)ID

……

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

with arcpy.da.SearchCursor(fc，["OID@"，"ORIG_FID"，"POINT_X"，"POINT_Y"，"Prev"，"Next"]，"OBJECTID="+str(sta)) as cursor： #獲取起點(diǎn)的信息

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點(diǎn)的信息

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點(diǎ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個(gè)夾角值的余弦值，并通過(guò)反三角函數(shù)求得第n個(gè)夾角的角度值。構(gòu)成該夾角的三個(gè)點(diǎn)分別為第n-1點(diǎn)，第n點(diǎn)，第n+1點(diǎn)，坐標(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)通過(guò)條件語(yǔ)句對(duì)第n個(gè)夾角余弦值與標(biāo)準(zhǔn)值的余弦值進(jìn)行比較，將構(gòu)成超限夾角的三個(gè)點(diǎn)狀數(shù)據(jù)導(dǎo)出到點(diǎn)問題圖層中，形成初步檢查結(jié)果。

3.4 生成問題線圖層

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

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

4.1 試驗(yàn)數(shù)據(jù)

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

4.2 試驗(yàn)過(guò)程

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

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

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

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

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

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

圖6 問題線圖層

4.3 試驗(yàn)分析

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

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

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

(3)效率更高，如果使用ArcGIS的折點(diǎn)夾角計(jì)算方法，需要耗費(fèi)的時(shí)間遠(yuǎn)比使用上述方法要長(zhǎng)。

5 總 結(jié)

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