基于ArcObject的扫描地形图坐标转换及自动化配准功能的实现

刘建民

(岳阳县林业局， 湖南 岳阳 414100)

基于ArcObject的扫描地形图坐标转换及自动化配准功能的实现

刘建民

(岳阳县林业局， 湖南 岳阳 414100)

采用VC++，基于ArcObject的组件二次开发模式，针对当前北京54坐标系地形图与西安80坐标系地形图共存的现状，通过对扫描地形图图幅四角坐标自动求算、不同坐标系间转换、扫描地形图配准等一系列自动化操作，极大简化了地形图电子化建库工作，实现了省级林业数据坐标系的统一。

地形图； 栅格配准； 坐标系； 坐标转换

近年湖南林业在数字化建设中，ArcGIS得到了较广泛的应用，也积累了众多林业电子地形图数据。由于历史原因，目前应用较广泛的是北京1954年坐标系地形图(简称北京54坐标系)，少量地方为西安1980年坐标系地形图(简称西安80坐标系)，也有些地方为两种坐标系地图混用的，有些地方甚至一直使用纸质地形图。

北京54坐标系是一种大地测量坐标系，其源于前苏联1942年普尔科夫坐标系，参考椭球为克拉索夫斯基椭球，由于当时的条件限制，北京54坐标系存在一些明显的缺点[1]，而西安80坐标系是通过对全国天文大地网整体平差而得，其采用多点定位原理建立，具有理论严密、定义明确、椭球参数精确、便于实际应用及理论研究、点位精度高等诸多特点[2]，使得西安80坐标近年来得到了全面的推广。随着林业信息化建设步伐的加快，特别是2014年湖南省林业厅下达《国家重点营造林工程小班数据库》建库指示精神，迫切需要一个自动化的地形图配准软件，完成统一坐标系的地形图基本数据建库工作。围绕这个目标，我们在VC++平台上利用ArcObject对象进行二次开发，实现标准分幅地形图四角坐标的求算、不同坐标系间坐标的转换、扫描地形图的配准与裁剪，全部操作基本自动化完成，为湖南省营造林数据库建库及林业数字化统一建设提供强大的技术支持。

1 坐标正反算

坐标正反算主要是指相同坐标系内进行地理坐标与高斯平面坐标之间的换算[3-4]。根据《国家基本比例尺地形图分幅和编号》，标准地形图一般是根据某一比例尺以一定的经差与纬差进行分幅的[5]，也就是说确定某张图幅，就同时确定了其图轮四角的经纬度坐标(地理坐标)。但在实际工作当中，我们一般用的是高斯平面坐标，因此使用时有必要进行换算操作。

1.1 大地坐标(B，L)转换为高斯平面坐标(x，y)-正算

x=X+Ntcos2B·l2×

y=NcosB·l×

式中：l为经度与中央子午线的经度差，t=tanB，η2=e′2cos2B，e′为地球第二偏心率，X为从赤道起算的子午线弧长[5-6]，计算公式为：

其中系数：

其中：e为椭球第一偏心率。

1.2 高斯平面坐标(x，y)转换为大地坐标(B，L)-反算

其中：Bf为底点纬度，下标“f”表示与Bf有关的量，底点纬度的计算公式是：

其中：

2 不同坐标系间坐标转换

从一个大地坐标系转换到另一个大地坐标系一般需要经过三个环节：大地坐标到地心坐标→地心坐标到地心坐标→地心坐标到大地坐标，同时涉及三种坐标类型：高斯平面直角坐标、大地坐标、空间直角坐标[7]。

2.1 大地坐标(B，L，H)转换为空间直角坐标(X，Y，Z)

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

2.2 空间直角坐标(X，Y，Z)转换为大地坐标(B，L，H)

已知(X，Y，Z)反算(B，L，H)时，采用迭代法实现，计算公式如下[7]：

2.3 赫尔默特(Helmert)转换求解七参数

坐标转换的中间环节，地心坐标到地心坐标，通常又称为7参数赫尔默特(Helmert)转换，将转换公式用7参数矩阵表示，即著名的布尔莎-活尔夫(Bursa-Wolf)公式[4，8-9]：

2.4 参数求算结果

参数求算最少要求3个点以上的公共点，以县域为单位的坐标系公共点的取得一般按本行政单位地域最外围四个方向分别取点。表1中，我们以湖南省某县为例分别选取其行政地域边界的5个公共点。求解七参数结果为：Δ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即以图幅号为基准进行标准图幅图轮四角原坐标取得，以及不同坐标系坐标转换的辅助界面。

图1 标准图幅图轮四角坐标转换软件界面Fig.1 Coordinate conversion software interface for standard images round corner

2.5 相关代码[10-11]

class CLjmCoordTrans7Param{

public：

CMatrix mMX；//7参数数组

public：

void Calc7Param(VCPnt3D& vcOld，

VCPnt3D& vcNew)；//计算7参数

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){

//创建A系数矩阵

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进行二次开发，所以可以直接利用ArcMap工具栏的“添加数据”工具即可。格栅图片文件添加之后，通过IMapPtr接口最终得到IRasterPtr用于栅格配准。

其中应用的接口与函数如下：

取得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]，并将屏幕坐标转换为像素坐标，并保存在数组变量中。

3.3 图轮四角地理坐标、平面直角坐标求算

一般依据《国家基本比例尺地形图分幅和编号》与图幅号可以直接求算本图幅西南角的地理坐标，然后依据标准图幅的经度差与纬度差，则可求得指定图幅图轮四角的地理坐标[5，13](见图1)。如果扫描地形图原有坐标系与要求配准的坐标系不同，则可先利用坐标转换模块将原有图轮的地理坐标转换为目标坐标系下的地理坐标。之后，再通过坐标正算代码将地理坐标(L，B)转换为高斯平面坐标(x，y)[4，11]。

3.4 投影参考的确定

投影参考的确定一般要用到ISpatialReferenceFactory2、ISpatialReference、IProjectedCoordinateSystem等接口[12]。创建投影参考的代码如下：

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 栅格图像配准

通过将鼠标在屏幕上拾取的四个图轮四角控制点按数据在X、Y轴方面上的大小关系可以实现对4个控制点的自动排序，然后求算获得的4个平面直角一一配对之后，利用IGeoReference接口、IRasterGeometryProc接口的Warp()与Register()函数即可完成栅格的配准[14]。

3.6 栅格图像的裁剪

由于扫描地形图在扫描之时保留了原纸质图的图轮四周外围空白部分，这将使得栅格配准之后会造成对相邻接图面的覆盖，因此必须切除多余部分才能将配准后的扫描地形图用于正式工作中。栅格图像的裁剪可以通过ArcObject的多个接口实现，如IExtractionOp接口[11]。关键代码如下：

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 栅格图层获取相关代码

获取第一个栅格图层代码如下：

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 结语

我们在VC++平台开发的《营造林制图软件》ArcForestry中集成了本文中的栅格自动化配准与坐标转换功能，完成一张地形图的配准工作只需轻松的用鼠标确定图轮四角即可，同时软件在ArcObject上开发，能够快捷便利的应用ArcGIS已有的界面资源，减少重复开发，加快了开发速度。

[1] 赵丹.基于VC++的常用大地坐标转换程序实现[J]. 铁道勘测与设计，2009(4)：11-17.

[2] 谢灵斌，韩广毅，丁孝兵.关于北京54与西安80坐标系相互转换的探讨[J]. 科技创新导报，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] 李岳.坐标转换系统的设计与实现[D]. 北京：中国地质大学，2010.

[8] 陈贻胜.坐标转换参数的求解方法及其应用[J]. 上海地质，2006(2)：53-56.

[9] 李潇，尹晖.基于最小二乘配置的三维空间坐标转换[J]. 测绘工程，2008，17(2)：16-17.

[10] 周长发.科学与工程数值算法[M]. 北京：清华大学出版社，2002.

[11] 牛丽娟.测量坐标转换模型研究与转换系统实现[D]. 西安：长安大学，2010.

[12] 牟乃夏，刘文宝.ArcGIS10地理信息系统教程[M]. 北京：测绘出版社，2012.

[13] 韩丽蓉.我国基本比例尺地形图分幅与编号的计算方法[J]. 青海大学学报(自然科学版)，2006，24(6)：79-82.

[14] 邱洪钢，张青莲，熊友谊.ArcGIS Engine地理信息系统开发[M]. 北京：人民邮电出版社，2013.

[15] 韩鹏，褚海峰.地理信息系统开发—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

国家林业科学数据平台(2005DKA32200-OG)。

刘建民(1973-)，男，湖南省岳阳县人，工程师，主要从事营造林、资源调查设计、林业信息化研究。

P 226.3

A

1003 — 5710(2015)03 — 0095 — 06

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

(文字编校：张 珉)