劉建民

(岳陽縣林業(yè)局， 湖南 岳陽 414100)

基于ArcObject的掃描地形圖坐標轉(zhuǎn)換及自動化配準功能的實現(xiàn)

劉建民

(岳陽縣林業(yè)局， 湖南 岳陽 414100)

采用VC++，基于ArcObject的組件二次開發(fā)模式，針對當前北京54坐標系地形圖與西安80坐標系地形圖共存的現(xiàn)狀，通過對掃描地形圖圖幅四角坐標自動求算、不同坐標系間轉(zhuǎn)換、掃描地形圖配準等一系列自動化操作，極大簡化了地形圖電子化建庫工作，實現(xiàn)了省級林業(yè)數(shù)據(jù)坐標系的統(tǒng)一。

地形圖； 柵格配準； 坐標系； 坐標轉(zhuǎn)換

近年湖南林業(yè)在數(shù)字化建設中，ArcGIS得到了較廣泛的應用，也積累了眾多林業(yè)電子地形圖數(shù)據(jù)。由于歷史原因，目前應用較廣泛的是北京1954年坐標系地形圖(簡稱北京54坐標系)，少量地方為西安1980年坐標系地形圖(簡稱西安80坐標系)，也有些地方為兩種坐標系地圖混用的，有些地方甚至一直使用紙質(zhì)地形圖。

北京54坐標系是一種大地測量坐標系，其源于前蘇聯(lián)1942年普爾科夫坐標系，參考橢球為克拉索夫斯基橢球，由于當時的條件限制，北京54坐標系存在一些明顯的缺點[1]，而西安80坐標系是通過對全國天文大地網(wǎng)整體平差而得，其采用多點定位原理建立，具有理論嚴密、定義明確、橢球參數(shù)精確、便于實際應用及理論研究、點位精度高等諸多特點[2]，使得西安80坐標近年來得到了全面的推廣。隨著林業(yè)信息化建設步伐的加快，特別是2014年湖南省林業(yè)廳下達《國家重點營造林工程小班數(shù)據(jù)庫》建庫指示精神，迫切需要一個自動化的地形圖配準軟件，完成統(tǒng)一坐標系的地形圖基本數(shù)據(jù)建庫工作。圍繞這個目標，我們在VC++平臺上利用ArcObject對象進行二次開發(fā)，實現(xiàn)標準分幅地形圖四角坐標的求算、不同坐標系間坐標的轉(zhuǎn)換、掃描地形圖的配準與裁剪，全部操作基本自動化完成，為湖南省營造林數(shù)據(jù)庫建庫及林業(yè)數(shù)字化統(tǒng)一建設提供強大的技術(shù)支持。

1 坐標正反算

坐標正反算主要是指相同坐標系內(nèi)進行地理坐標與高斯平面坐標之間的換算[3-4]。根據(jù)《國家基本比例尺地形圖分幅和編號》，標準地形圖一般是根據(jù)某一比例尺以一定的經(jīng)差與緯差進行分幅的[5]，也就是說確定某張圖幅，就同時確定了其圖輪四角的經(jīng)緯度坐標(地理坐標)。但在實際工作當中，我們一般用的是高斯平面坐標，因此使用時有必要進行換算操作。

1.1 大地坐標(B，L)轉(zhuǎn)換為高斯平面坐標(x，y)-正算

x=X+Ntcos2B·l2×

y=NcosB·l×

式中：l為經(jīng)度與中央子午線的經(jīng)度差，t=tanB，η2=e′2cos2B，e′為地球第二偏心率，X為從赤道起算的子午線弧長[5-6]，計算公式為：

其中系數(shù)：

其中：e為橢球第一偏心率。

1.2 高斯平面坐標(x，y)轉(zhuǎn)換為大地坐標(B，L)-反算

其中：Bf為底點緯度，下標“f”表示與Bf有關(guān)的量，底點緯度的計算公式是：

其中：

2 不同坐標系間坐標轉(zhuǎn)換

從一個大地坐標系轉(zhuǎn)換到另一個大地坐標系一般需要經(jīng)過三個環(huán)節(jié)：大地坐標到地心坐標→地心坐標到地心坐標→地心坐標到大地坐標，同時涉及三種坐標類型：高斯平面直角坐標、大地坐標、空間直角坐標[7]。

2.1 大地坐標(B，L，H)轉(zhuǎn)換為空間直角坐標(X，Y，Z)

式中：e2為第一偏心率平方，N為卯酉圈曲率半徑。

2.2 空間直角坐標(X，Y，Z)轉(zhuǎn)換為大地坐標(B，L，H)

已知(X，Y，Z)反算(B，L，H)時，采用迭代法實現(xiàn)，計算公式如下[7]：

2.3 赫爾默特(Helmert)轉(zhuǎn)換求解七參數(shù)

坐標轉(zhuǎn)換的中間環(huán)節(jié)，地心坐標到地心坐標，通常又稱為7參數(shù)赫爾默特(Helmert)轉(zhuǎn)換，將轉(zhuǎn)換公式用7參數(shù)矩陣表示，即著名的布爾莎-活爾夫(Bursa-Wolf)公式[4，8-9]：

2.4 參數(shù)求算結(jié)果

參數(shù)求算最少要求3個點以上的公共點，以縣域為單位的坐標系公共點的取得一般按本行政單位地域最外圍四個方向分別取點。表1中，我們以湖南省某縣為例分別選取其行政地域邊界的5個公共點。求解七參數(shù)結(jié)果為：Δx=-310.361 305，Δy=784.770 884，Δz=562.618 404，m=-0.000175，Rx=0.000 169，Ry=-0.000 368，Rz=-0.000 218。

表1 坐標系A,B中的公共點Tab.1 CommonpointinthecoordinatesystemofAandB點號坐標系AXYZ1462823.03206764.0374.02451240.03238787.0908.23401000.03210000.039.04403000.03210000.043.05381000.03265000.025.3點號坐標系BXYZ1462771.63206795.2373.22451184.83238728.2908.23400940.93209949.039.04402939.93209946.943.05380941.93264946.525.3

圖1即以圖幅號為基準進行標準圖幅圖輪四角原坐標取得，以及不同坐標系坐標轉(zhuǎn)換的輔助界面。

圖1 標準圖幅圖輪四角坐標轉(zhuǎn)換軟件界面Fig.1 Coordinate conversion software interface for standard images round corner

2.5 相關(guān)代碼[10-11]

class CLjmCoordTrans7Param{

public：

CMatrix mMX；//7參數(shù)數(shù)組

public：

void Calc7Param(VCPnt3D& vcOld，

VCPnt3D& vcNew)；//計算7參數(shù)

void conv3dOld2New(VCPnt3D& vcOld，

VCPnt3D& vcNew)；

}；

void CLjmCoordTrans7Param：：Calc7Param(

VCPnt3D& vcOld，VCPnt3D& vcNew){

int iCnt=vcOld.size()；

CMatrix txA；

CMatrix txL；

CMatrix txAT(7，3*iCnt)；

CMatrix txATA(7，7)；

CMatrix txInvATA；

CMatrix txATL(7，1)；

txA=GetMxA(vcOld)；

txL=GetMxL(vcOld，vcNew)；

txAT=txA.Transpose()；

txATA=txAT*txA；

txATL=txAT*txL；

txInvATA=txATA；

txInvATA.InvertGaussJordan()；

mMX=txInvATA*txATL；

}

void CLjmCoordTrans7Param：：conv3dOld2New(

VCPnt3D& vcOld，VCPnt3D& vcNew){

//創(chuàng)建A系數(shù)矩陣

CMatrix mxA(iCnt*3，7)；

CMatrix mxAi(3，7)；

CMatrix mxS(iCnt*3，1)；

CLjmPoint3D p3d；

for(i=0；i

p3d=vcOld.at(i)；

mxAi=Trans3DToAi(p3d)；

PutMatrA(mxA，mxAi，i)；

mxS.SetElement(i*3+0，0，p3d.mX)；

mxS.SetElement(i*3+1，0，p3d.mY)；

mxS.SetElement(i*3+2，0，p3d.mZ)；

}

CMatrix mxT；

mxT=(mxA*mMX)+mxS；

iCnt=mxT.GetNumRows()/3；

vcNew.clear()；

vcNew.erase(vcNew.begin()，vcNew.end())；

for(i=0；i

p3d.Set(mxT.GetElement(i*3，0)，

mxT.GetElement(i*3+1，0)，

mxT.GetElement(i*3+2，0))；

vcNew.push_back(p3d)；

}

}

3 柵格配準

3.1 加載格柵

由于本軟件是在ArcGIS基礎上利用ArcObject進行二次開發(fā)，所以可以直接利用ArcMap工具欄的“添加數(shù)據(jù)”工具即可。格柵圖片文件添加之后，通過IMapPtr接口最終得到IRasterPtr用于柵格配準。

其中應用的接口與函數(shù)如下：

取得IEnumLayer對象：

IMap：：get_Layers(NULL，VARIANT_FALSE，IEnumLayer** ipEL)

取得第一個柵格圖層：

IEnumLayer：：Next(ILayer** ipL)

取得IRaster 對象：

IRasterLayer：：get_Raster(IRaster** ipR)；

3.2 拾取控制點

標準圖幅圖輪四角控制點格柵像素坐標的拾取一般通過ArcObject的ITool接口從OnMouseDown事件上獲得屏幕坐標[12]，并將屏幕坐標轉(zhuǎn)換為像素坐標，并保存在數(shù)組變量中。

3.3 圖輪四角地理坐標、平面直角坐標求算

一般依據(jù)《國家基本比例尺地形圖分幅和編號》與圖幅號可以直接求算本圖幅西南角的地理坐標，然后依據(jù)標準圖幅的經(jīng)度差與緯度差，則可求得指定圖幅圖輪四角的地理坐標[5，13](見圖1)。如果掃描地形圖原有坐標系與要求配準的坐標系不同，則可先利用坐標轉(zhuǎn)換模塊將原有圖輪的地理坐標轉(zhuǎn)換為目標坐標系下的地理坐標。之后，再通過坐標正算代碼將地理坐標(L，B)轉(zhuǎn)換為高斯平面坐標(x，y)[4，11]。

3.4 投影參考的確定

投影參考的確定一般要用到ISpatialReferenceFactory2、ISpatialReference、IProjectedCoordinateSystem等接口[12]。創(chuàng)建投影參考的代碼如下：

ISpatialReferencePtr ipSRef；

long lProjCS；

lProjCS=coo.TransSRProjCS4TypeCode_JW(bHaveDiaHao)；

ISpatialReferenceFactory2Ptr ipSRFac；

ipSRFac.CreateInstance(CLSID_SpatialReferenceEnvironment)；

IProjectedCoordinateSystemPtr ipProjCS；

ipSRFac->CreateProjectedCoordinateSystem(lProjCS，& ip ProjCS)；

ipSRef=ipProjCS；

其中投影代號計算如下(以3度帶為例)：

long lProjCS =0；

if(eERef==fctERef_BJ54){

lProjCS =iDaiHao -25 +2422；

}else if(eERef==fctERef_XIAN80){

lProjCS =iDaiHao -25 +2370；

}

3.5 柵格圖像配準

通過將鼠標在屏幕上拾取的四個圖輪四角控制點按數(shù)據(jù)在X、Y軸方面上的大小關(guān)系可以實現(xiàn)對4個控制點的自動排序，然后求算獲得的4個平面直角一一配對之后，利用IGeoReference接口、IRasterGeometryProc接口的Warp()與Register()函數(shù)即可完成柵格的配準[14]。

3.6 柵格圖像的裁剪

由于掃描地形圖在掃描之時保留了原紙質(zhì)圖的圖輪四周外圍空白部分，這將使得柵格配準之后會造成對相鄰接圖面的覆蓋，因此必須切除多余部分才能將配準后的掃描地形圖用于正式工作中。柵格圖像的裁剪可以通過ArcObject的多個接口實現(xiàn)，如IExtractionOp接口[11]。關(guān)鍵代碼如下：

IExtractionOpPtr ipEXOP(CLSID_RasterExtractionOp)；

IRasterPtr ipRsrc=GetRaster()；

IGeoDatasetPtr ipGDSsrc(ipRsrc)；

VARIANT vVal=GetCellSizeVarX()；

IRasterAnalysisEnvironmentPtr ipRAEn =ipEXOP；

ipRAEn->SetCellSize(esriRasterEnvValue，&vVal)；

IEnvelopePtr ipEnv；

ipPgnClip->get_Envelope(&ipEnv)；

VARIANT vExtProvider =_variant_t(ipEnv，true)；

ipRAEn->SetExtent(esriRasterEnvValue，

&vExtProvider，NULL)；

IGeoDatasetPtr ipGDSaim；

IRasterLayerPtr ipRLaim；

IRasterDatasetPtr ipRDSaim；

IRasterPtr ipRaim；

ipEXOP->Polygon(ipGDSsrc，ipPgnClip，

VARIANT_TRUE，&ipGDSaim)；

也可以應用IConversionOp等接口進行裁剪操作。

3.7 柵格圖層獲取相關(guān)代碼

獲取第一個柵格圖層代碼如下：

iLayerPtr ipL=NULL；

IRasterLayerPtr ipRL=NULL；

IEnumLayerPtr ipEL；

ipMap->get_Layers(NULL，VARIANT_FALSE，& ip EL)；

if(ipEL){

ipEL->Next(&ipL)；

while(ipL){

ipRL=ipL；

if(ipRL)break；

ipEL->Next(&ipL)；

}

}

4 程序流程圖[15]

5 結(jié)語

我們在VC++平臺開發(fā)的《營造林制圖軟件》ArcForestry中集成了本文中的柵格自動化配準與坐標轉(zhuǎn)換功能，完成一張地形圖的配準工作只需輕松的用鼠標確定圖輪四角即可，同時軟件在ArcObject上開發(fā)，能夠快捷便利的應用ArcGIS已有的界面資源，減少重復開發(fā)，加快了開發(fā)速度。

[1] 趙丹.基于VC++的常用大地坐標轉(zhuǎn)換程序?qū)崿F(xiàn)[J]. 鐵道勘測與設計，2009(4)：11-17.

[2] 謝靈斌，韓廣毅，丁孝兵.關(guān)于北京54與西安80坐標系相互轉(zhuǎn)換的探討[J]. 科技創(chuàng)新導報，2010(16)：16-88.

[3] 董金壯，趙淵新，刁芹元.高斯投影變換程序的編寫[J]. 新疆有色金屬，2003(增刊)：8-9.

[4] 孔祥元，郭際明，劉宗泉.大地測量學基礎[M]. 武漢：武漢大學出版社，2005.

[5] GB/T 13989-92，國家基本比例尺地形圖分幅和編號[S].

[6] 劉正才.子午線弧長公式的簡化及通用高斯投影計算程序介紹[J]. 測繪工程，2001，10(1)：55-57.

[7] 李岳.坐標轉(zhuǎn)換系統(tǒng)的設計與實現(xiàn)[D]. 北京：中國地質(zhì)大學，2010.

[8] 陳貽勝.坐標轉(zhuǎn)換參數(shù)的求解方法及其應用[J]. 上海地質(zhì)，2006(2)：53-56.

[9] 李瀟，尹暉.基于最小二乘配置的三維空間坐標轉(zhuǎn)換[J]. 測繪工程，2008，17(2)：16-17.

[10] 周長發(fā).科學與工程數(shù)值算法[M]. 北京：清華大學出版社，2002.

[11] 牛麗娟.測量坐標轉(zhuǎn)換模型研究與轉(zhuǎn)換系統(tǒng)實現(xiàn)[D]. 西安：長安大學，2010.

[12] 牟乃夏，劉文寶.ArcGIS10地理信息系統(tǒng)教程[M]. 北京：測繪出版社，2012.

[13] 韓麗蓉.我國基本比例尺地形圖分幅與編號的計算方法[J]. 青海大學學報(自然科學版)，2006，24(6)：79-82.

[14] 邱洪鋼，張青蓮，熊友誼.ArcGIS Engine地理信息系統(tǒng)開發(fā)[M]. 北京：人民郵電出版社，2013.

[15] 韓鵬，褚海峰.地理信息系統(tǒng)開發(fā)—ArcObjects方法[M]. 武漢：武漢大學出版社，2005.

RealizationofscannedtopographicmapcoordinateconversionandautomaticregistrationfunctionbasedonArcObject

LIU Jianmin

(Forestry Bureau of Yueyang County， Yueyang 414100， China)

Used VC++IDE，based on ArcObject of component II times development mode，for current Beijing 54 and Xi’an 80 coordinates topographic maps coexistence of status，through a series of automatic operation on automatically calculation of four corners coordinates based on scanned topographic maps，and conversion between different coordinates systems，and registration of scanned topographic maps，the work of topographic maps electronic database was greatly simplified，the unity of the provincial forestry data coordinate system was achieved.

topographic map； raster registration； coordinate system； coordinate transformation

2015-03-13

國家林業(yè)科學數(shù)據(jù)平臺(2005DKA32200-OG)。

劉建民(1973-)，男，湖南省岳陽縣人，工程師，主要從事營造林、資源調(diào)查設計、林業(yè)信息化研究。

P 226.3

A

1003 — 5710(2015)03 — 0095 — 06

10. 3969/j. issn. 1003 — 5710. 2015. 03. 021

(文字編校：張 珉)