閭海洋,湯 玲,王 剛
(江蘇省地質(zhì)測繪院,江蘇 南京 211102)
國務(wù)院于2017年10月16日下發(fā)《國務(wù)院關(guān)于開展第三次全國土地調(diào)查的通知》(國發(fā)〔2017〕48號),決定自2017年起開展第三次全國土地調(diào)查(以下簡稱“三調(diào)”)。為統(tǒng)一三調(diào)技術(shù)標(biāo)準(zhǔn),方便各個(gè)縣區(qū)的數(shù)據(jù)庫能夠融合到一起,自然資源部于2019年1月28日發(fā)布了《第三次全國國土調(diào)查技術(shù)規(guī)程》(TD/T 1055-2019,以下簡稱《三調(diào)技術(shù)規(guī)程》)。三調(diào)以縣(區(qū))為單位進(jìn)行調(diào)查,每個(gè)縣都有一個(gè)控制面積,每個(gè)縣內(nèi)所有圖斑的面積和等于縣控制面積,所有縣的面積和等于國家的控制面積。通過這種方式進(jìn)行控制,能夠保證全國的面積是一個(gè)固定值,不會因?yàn)橛?jì)算方式的不同而產(chǎn)生不同的總面積。三調(diào)圖斑沒有采用高斯投影面上圖形的幾何面積,而是使用橢球面積作為圖斑的面積[1]。計(jì)算圖斑橢球面積這一功能相對復(fù)雜,在主流的GIS建庫軟件中都沒有現(xiàn)成的工具,因此需要使用GIS軟件提供的開發(fā)接口進(jìn)行二次開發(fā)。
《三調(diào)技術(shù)規(guī)程》中對圖斑的橢球面積計(jì)算方式做出了明確的規(guī)定:對任意一個(gè)封閉區(qū)域(圖斑)進(jìn)行分割,將封閉區(qū)域分割成多個(gè)梯形圖塊,先對每個(gè)梯形圖塊進(jìn)行計(jì)算,然后再對所有梯形的面積求和從而得到封閉區(qū)域的面積。對于橢球面上的任意一個(gè)梯形圖塊的應(yīng)采用公式(1)進(jìn)行計(jì)算:
(1)
其中,a表示橢球的長半軸,單位為m;b表示橢球的短半軸,單位為m;ΔL表示梯形圖塊的經(jīng)度差值,單位為弧度;BΔ表示梯形圖塊的緯度差值,單位為弧度;Bm表示梯形圖塊的南北緯度的平均值,單位為弧度。
《三調(diào)技術(shù)規(guī)程》規(guī)定平面坐標(biāo)使用2000國家大地坐標(biāo)系,因此參數(shù)α、b、e2、A、B、C、D、E都為常數(shù),在計(jì)算圖斑面積之前可以固化在代碼中,具體數(shù)值如表1所示。
表1 基本常數(shù)表
圖斑是土地實(shí)際的利用狀況在數(shù)據(jù)庫的表現(xiàn),都是不規(guī)則的多邊形,因此需要對圖斑進(jìn)行分割處理,將不規(guī)則的圖斑分割成若干個(gè)規(guī)則的梯形圖塊進(jìn)行計(jì)算,分割方式如圖1所示。
圖1 任意圖斑梯形分割示意圖
因此圖斑ABCDE的總面積公式可以表示為公式(2):
SABCDE=SABB1A1+SBCC1B1+SCDD1C11+SDEE1D1+SEAA1E1
(2)
式中的每個(gè)面積子項(xiàng)都可以通過梯形圖塊面積計(jì)算公式得到。計(jì)算梯形圖塊面積時(shí),使用的坐標(biāo)是地理坐標(biāo),而實(shí)際生產(chǎn)過程中使用的坐標(biāo)是高斯平面坐標(biāo),在計(jì)算時(shí)還需要通過高斯投影反算公式將高斯平面坐標(biāo)反算為地理坐標(biāo)。
在計(jì)算圖斑的橢球面積過程中,為提高面積計(jì)算精度,如果圖斑節(jié)點(diǎn)上任意相鄰節(jié)點(diǎn)之間的線段長度大于70 m時(shí),需要進(jìn)行插值計(jì)算。內(nèi)插點(diǎn)僅作為橢球面積計(jì)算使用,不需要對線段進(jìn)行分割。進(jìn)行內(nèi)插時(shí),內(nèi)插點(diǎn)應(yīng)分布均勻,插值計(jì)算流程如圖2所示。
圖2 插值計(jì)算流程
通過面積計(jì)算公式計(jì)算出每一個(gè)圖斑的橢球面積,但由于計(jì)算過程中會產(chǎn)生一些誤差,最終導(dǎo)致面積匯總時(shí)總面積會出現(xiàn)不確定性。為解決這一問題,三調(diào)采用了控制面積制度。國土調(diào)查數(shù)據(jù)庫中的縣級及縣級以上行政區(qū)域界線采用全國陸地行政區(qū)域勘界成果確定的界線,鄉(xiāng)鎮(zhèn)級行政區(qū)域界線采用各縣(市、區(qū))最新確定的界線。為保證全國各個(gè)縣(市、區(qū))的國土調(diào)查數(shù)據(jù)庫拼到一起后的總面積的確定性,在計(jì)算面積時(shí)實(shí)行了控制面積制度。在開展三調(diào)建庫時(shí),國家統(tǒng)一下發(fā)了各個(gè)縣(市、區(qū))的控制面積,在一個(gè)縣(市、區(qū))范圍內(nèi),所有圖斑的總面積應(yīng)等于下發(fā)控制面積。面積計(jì)算以平方米為單位,保留2位有效小數(shù)。
在三調(diào)中,圖斑號以村為單位進(jìn)行編號,控制面積實(shí)行三級控制,即縣控制各個(gè)鄉(xiāng)鎮(zhèn)的面積,鄉(xiāng)鎮(zhèn)控制各個(gè)村(村級調(diào)查區(qū)),村再控制內(nèi)部的各個(gè)圖斑。面積平差的總體流程為:① 計(jì)算各個(gè)鄉(xiāng)鎮(zhèn)的橢球面積;② 匯總各個(gè)鄉(xiāng)鎮(zhèn)面積,得到匯總面積與縣級控制面積的差值,再將差值以鄉(xiāng)鎮(zhèn)面積大小為權(quán)分配到各個(gè)鄉(xiāng)鎮(zhèn)中;③ 計(jì)算各個(gè)村(村級調(diào)查區(qū))的橢球面積;④ 匯總各個(gè)村的面積,得到匯總面積與鎮(zhèn)控制面積的差值,再將差值以村面積大小為權(quán)分配到各個(gè)村中;⑤ 計(jì)算各個(gè)圖斑的橢球面積;⑥ 以村為單位進(jìn)行圖斑橢球面積匯總,得到匯總面積與村控制面積的差值,再將差值以圖斑面積大小為權(quán)分配到各個(gè)圖斑中。
橢球面積計(jì)算工具分為3個(gè)模塊來實(shí)現(xiàn),分別為數(shù)據(jù)管理模塊、橢球面積計(jì)算模塊和面積平差模塊。3個(gè)模塊相互調(diào)用,共同實(shí)現(xiàn)圖斑橢球面積計(jì)算的全部功能。
三調(diào)數(shù)據(jù)庫是一個(gè)空間數(shù)據(jù)庫中,行政區(qū)(鎮(zhèn))、村級調(diào)查區(qū)與地類圖斑分別保存于XZQ、CJDCQ、DLTB層中,數(shù)據(jù)結(jié)構(gòu)如表2-表4所示[2]。
表2 XZQ屬性結(jié)構(gòu)表(部分)
表3 CJDCQ屬性結(jié)構(gòu)表(部分)
表4 DLTB屬性結(jié)構(gòu)表(部分)
此模塊主要功能有:行政區(qū)、村級調(diào)查區(qū)與地類圖斑數(shù)據(jù)的加載、讀取、賦值,要素類屬性結(jié)構(gòu)檢查,要素的遍歷、空間參考獲取等操作。使用IWorkspaceFactory.Open方法將三調(diào)數(shù)據(jù)庫加載到內(nèi)存,得到Workspace對象,Workspace對象同時(shí)實(shí)現(xiàn)了IFeatureWorkspace。使用IFeatureWorkspace. OpenFeatureClass方法代碼行政區(qū)或地類圖斑的名稱,得到FeatureClass對象。然后通過IFeatureClass.Search方法獲取IFeatureCursor對象后,就可以利用IFeatureCursor.NextFeature方法實(shí)現(xiàn)對行政區(qū)、地類圖斑層中各個(gè)要素的遍歷[3]。
行政區(qū)、村級調(diào)查區(qū)與地類圖斑在空間數(shù)據(jù)庫中都是以多邊形(Polygon)的形式存在,可以使用統(tǒng)一的方法計(jì)算其橢球面積。在GeoDatabase數(shù)據(jù)模型中,多邊形是由一個(gè)或多個(gè)環(huán)(Ring)組成的。環(huán)可以分為外環(huán)(ExteriorRing)和內(nèi)環(huán)(InteriorRing),一個(gè)外環(huán)內(nèi)可以有0個(gè)或多個(gè)內(nèi)環(huán),內(nèi)環(huán)不能獨(dú)立存在。外環(huán)的面積為正數(shù),內(nèi)環(huán)的面積為負(fù)數(shù)。對于沒有島洞的多邊形而言,其只有一個(gè)環(huán)[4-5]。
橢球面積計(jì)算模塊主要功能是計(jì)算環(huán)的橢球面積及多邊形內(nèi)部環(huán)的遍歷。計(jì)算環(huán)的橢球面積是此模塊中的重點(diǎn),也是整個(gè)工具中最難實(shí)現(xiàn)的部分。
對于一個(gè)環(huán)而言,其主要接口為IRing,但環(huán)對象也實(shí)現(xiàn)了IPointCollection接口,通過對環(huán)中每個(gè)節(jié)點(diǎn)的遍歷,每兩個(gè)連續(xù)點(diǎn)組成一個(gè)線段,實(shí)現(xiàn)橢球面積的計(jì)算,主要實(shí)現(xiàn)方法為:① 初始化節(jié)點(diǎn)索引N=0 ;② 通過IPointCollection.Point方法,代入索引號得到邊 Edeg(N,N+1);③ 計(jì)算Edeg(N,N+1)的長度,如果小于或等于70 m,則直接使用,如果大于70 m,則進(jìn)行內(nèi)插,形成一個(gè)或多個(gè)子邊;④ 對于每個(gè)子邊而言,將首尾結(jié)點(diǎn)通過IPoint.Project方法計(jì)算點(diǎn)的經(jīng)緯度坐標(biāo);⑤ 通過梯形面積公式(1)計(jì)算每個(gè)梯形的面積;⑥ 使用N=N+1,回到步驟2再進(jìn)行計(jì)算,直到計(jì)算到最后一個(gè)點(diǎn);⑦ 使用公式(2)中提供的方法,匯總所有梯形的面積,得到環(huán)的面積。
XZQ層的橢球面積保存于JSMJ字段中,CJDCQ層的橢球面積保存于JSMJ字段中,DLTB層的橢球面積保存于TBMJ字段中。
面積平差關(guān)系到整個(gè)數(shù)據(jù)庫的面積匯總值能否與國家下發(fā)的控制面積保持一致。受計(jì)算策略、計(jì)算精度及小位取位的影響,直接計(jì)算得到的面積很難與控制面積保持一致,因此得到圖斑的橢球面積后,必須要進(jìn)行面積平差。進(jìn)行平差時(shí),分配多出面積的基本策略就以橢球面積值為權(quán)進(jìn)行分配,橢球面積值越大,分配到的誤差值就越大。關(guān)鍵算法如下:
(1)匯總XZQ層的JSMJ字段,得到總控制面積與匯總面積的差值,將差值以JSMJ值為權(quán)進(jìn)行分配,平差后的面積保存于DCMJ字段。
(2)一個(gè)XZQ內(nèi)部有一個(gè)或多個(gè)CJDCQ,以XZQ為單位,分組匯總CJDCQ的JSMJ字段,得到每組的匯總值與相應(yīng)XZQ的DCMJ值之間的差值,再將差值以CJDCQ的JSMJ值為權(quán)進(jìn)行分配,平差后的面積值保存于CJDCQ的DCMJ字段。
(3)一個(gè)CJDCQ內(nèi)部有一個(gè)或多個(gè)DLTB,以CJDCQ為單位,分組匯總DLTB的TBMJ,得到每組的匯總值與相應(yīng)CJDCQ的DCMJ值之間的差值,再將差值以DLTB的TBMJ值為權(quán)進(jìn)行分配,平差后的面積值保存于DLTB的TBMJ字段。
在進(jìn)行面積平差時(shí),受小數(shù)取位的影響,可能一次分配不能達(dá)到匯總值與控制值完全相同,還需要進(jìn)行第二次分配,再次分配還如果還有差值,那么就將這個(gè)差值直接分配給面積最大的那個(gè)圖形。
橢球面積計(jì)算工具的開發(fā)基于ArcGIS Add-in框架,IDE選擇Visual Studio2013,使用VisualBasic.NET作為編程語言實(shí)現(xiàn)以上的全部功能。
所有代碼編寫、調(diào)試完成后,即可對此項(xiàng)目進(jìn)行編譯。編譯通過后,會生成一個(gè)后綴為esriAddIn的安裝包文件。將此文件拷到具有相應(yīng)ArcGIS版本的電腦上,雙擊此文件后在彈出的窗口中點(diǎn)擊“安裝”完成工具部署。打開ArcMap軟件,在工具欄中將此工具條顯示出來即可使用。圖3為橢球面積計(jì)算工具的操作界面。
第三次國土調(diào)查已經(jīng)在全國范圍內(nèi)開展,橢球面積計(jì)算工具已經(jīng)應(yīng)用到了我院承接了多個(gè)區(qū)縣的三調(diào)項(xiàng)目中,使用此工具能夠快速對縣級數(shù)據(jù)庫進(jìn)行面積計(jì)算,計(jì)算結(jié)果能夠通過國家質(zhì)檢軟件的檢查。
圖3 橢球面積計(jì)算
本文通過橢球面積計(jì)算方法的分析,并以ArcMap AddIn為開發(fā)框架,開發(fā)了專用的橢球面積計(jì)算工具,能夠在ArcMap中快速進(jìn)行橢球面積計(jì)算,方便內(nèi)業(yè)人員的建庫工作。文中關(guān)于ArcMap Add-in開發(fā)的案例對開發(fā)其他基于ArcMap的應(yīng)用也有一定的借鑒意義。