基于ArcGIS Engine的矢柵基礎地理數據空間轉換方法研究

2011-11-15 03:29:52劉佰瓊戴相喜

測繪通報 2011年12期

劉佰瓊，周衛，戴相喜

(1.南京師范大學教育部虛擬地理環境重點實驗室，江蘇南京210046;2.南京市測繪勘察研究院有限公司，江蘇南京210005)

劉佰瓊1，周衛1，戴相喜2

(1.南京師范大學教育部虛擬地理環境重點實驗室，江蘇南京210046;2.南京市測繪勘察研究院有限公司，江蘇南京210005)

隨著2000國家大地坐標系的正式啟用，實現基礎地理數據不同坐標系間的空間轉換具有重要意義和實用價值。基于Esri二次開發組件ArcGIS Engine所提供的空間轉換和影像配準接口，討論ArcGIS平臺坐標轉換和投影變換方法，采用C#語言設計開發空間轉換系統，并進行轉換誤差分析，驗證該轉換方法的實用性和可靠性。

空間轉換;投影變換;ArcGIS Engine

一、引言

經國務院批準，我國自2008年7月1日起，正式啟用2000國家大地坐標系(簡稱CGCS2000)［1-2］。然而我國現有的測繪成果大多是基于1954北京坐標系、1980西安坐標系或地方獨立坐標系獲得的，因此在今后相當長的時期內，面臨著大量的基礎地理數據空間轉換問題。

目前許多學者都對空間轉換模型進行了研究，探討了公共點的選取原則和模型適用范圍［3-4］。本文針對GIS中矢量數據和柵格數據的結構特點，討論了ArcGIS軟件平臺坐標系統和地理參考投影變換的相關接口，并以實例介紹了在 Visual Studio 2005集成開發環境下使用Esri二次開發組件Arc-GIS Engine實現矢柵基礎地理數據空間轉換的關鍵技術方法。

二、矢量數據空間轉換

1.轉換模型

二維的平面坐標轉換模型有多項式逼近、正形轉換、相似變換、仿射變換等，其中相似變換四參數模型具有較高的嚴密性，適用于區域范圍小、轉換區域在兩個坐標系中的幾何相似程度好的情況［5-6］，本文采用平面四參數模型進行矢量數據空間轉換。

2.轉換方法與流程

矢量數據是基礎地理數據最主要的數據類型，能準確表示地理空間實體的位置、長度、面積等特性。由于矢量數據的所有空間幾何實體都是由坐標點(串)構成的，所以逐點進行空間轉換是最嚴密的轉換方式。但是逐點轉換需要對所有圖形要素重新生成，要求讀取并保存不同數據類型的符號，并且必須保證轉換前后符號的一致性，增大了系統開銷，計算量大，運行效率低。為了在提高轉換效率的同時保證轉換成果的可靠性，本文采用逐要素轉換的方法，對矢量數據的每個圖元要素依次進行轉換。矢量數據的空間轉換流程如圖1所示，其中包括數據輸入、數據處理和數據輸出3大模塊。

圖1 矢量數據空間轉換流程圖

3.轉換功能的實現

在二次開發組件 ArcGIS Engine中，提供了ITransform2D接口實現圖形要素的平移、旋轉和縮放，還提供了強大的坐標轉換方法(transform)。ITransformation接口則提供了多種轉換方法，包括莫洛金斯基轉換(Abridged Molodensky transformation)、二次多項式轉換(Affine transformation)等。矢量數據空間轉換的主要代碼與說明如下:

(1)根據控制點坐標擬合四參數

輸入已知的n組控制點坐標，根據這些控制點坐標采用最小二乘擬合的方法計算轉換參數。

IPoint［］fromPts=new PointClass［count］;

IPoint［］toPts=new PointClass［count］;

∥定義并初始化點數組，其中count為控制點個數

IAffineTransformation2D3GEN pAffineTrans=new Affine-Transformation2DClass();

∥定義并實例化相似變換接口對象

pAffineTrans. DefineConformalFromControlPoints (ref fromPts，ref toPts);

∥計算轉換參數，其中Dx、Dy為坐標平移量，Dxscale、Dyscale為縮放尺度，Drotate為旋轉角度

Dx=pAffineTrans.XTranslation;

Dy=pAffineTrans.YTranslation;

Dxscale=pAffineTrans.XScale;

Dyscale=pAffineTrans.YScale;

Drotate=pAffineTrans.Rotation;

(2)圖元要素的空間轉換

ArcGIS中的圖元要素分為點、線、面、注記4種類型，其中點要素由坐標點對構成，線要素由一系列坐標串構成，面要素則由一系列首尾相連的坐標串構成。點、線、面3類要素都是一般的圖形要素類，只需對其進行平移、旋轉、縮放就可實現空間轉換。

IFeatureCursor pFeaCur=pFeaCls.Update(null，false);

IFeature pFeature=pFeaCur.NextFeature();

IPoint pPoint=new PointClass();∥坐標基點

pPoint.X=Dx;pPoint.Y=Dy;

ITransform2D pTransform2D=pFeature.Shape;

∥要素按照四參數進行平移、旋轉和縮放

pTransform2D.Move(Dx，Dy);

pTransform2D.Rotate(pPoint，Drotate);

pTransform2D.Scale(pPoint，Dxscale，Dyscale);

pFeature.Shape=pTransform2D as IGeometry;

∥將轉換后的要素更新保存

pFeatureCursor.UpdateFeature(pFeature);

注記要素包含文字的空間位置、字體等注記信息，要求空間轉換之后保持字體朝向不變，因此可采用中心點位平移的方法進行注記的空間轉換。即首先計算出注記中心點經四參數模型轉換后的坐標;再根據其轉換前后的坐標值計算出坐標偏移量;最后將注記的文字和圖形均按照中心點偏移量進行空間轉換。

三、柵格數據空間轉換

1.轉換方法與流程

一般而言，要將柵格數據的空間信息轉換到其他坐標系中，可通過幾何校正來實現［7］。但是幾何校正需要人工選取同名控制點，轉換效率比較低。本文探討了在已知兩個坐標系之間轉換參數的情況下，自動實現柵格數據空間轉換的方法。

柵格數據空間轉換流程如圖2所示，首先獲得柵格數據的四至坐標，并將其作為空間轉換的控制點。然后通過已知的轉換參數計算出其在目標坐標系下的新坐標，再在此基礎上對柵格像元進行重采樣，實現空間轉換。該方法無需人工選取控制點，根據轉換參數就可以方便地實現柵格數據不同坐標系之間的空間轉換，提高了柵格數據空間轉換的自動化程度，尤其適用于批量轉換。

圖2 柵格數據坐標轉換流程圖

2.轉換功能的實現

柵格數據空間轉換可通過IRaster GeometryProc接口實現，該接口提供了柵格數據的剪切(clip)、翻轉(filp)、合并(merge)、映像(mirror)及鑲嵌(mosaic)等操作。通過轉換參數對柵格數據進行空間轉換可以使用該接口的Wrap方法，該方法會略微改變影像圖的色彩值。柵格數據空間轉換的主要代碼與說明如下:

1)獲取柵格數據源坐標系的四至坐標

IEnvelope pEnv=pRasterLayer.AreaOfInterest;

∥根據最小外接矩形獲得柵格數據四至坐標，并將其存入fromPtsCol點數組，作為源控制點坐標

IPointCollection fromPtsCol=new MultipointClass();

IPoint frPt1=new PointClass();

frPt1.X=pEnv.XMin;frPt1.Y=pEnv.YMin;

fromPtsCol.AddPoint(frPt1，ref obj，ref obj);

……

2)計算柵格數據目標坐標系的四至坐標

∥將四至坐標根據四參數進行坐標轉換，并將轉換結果存入toPtsCol點數組，作為目標控制點坐標

IClone pClone=frPt1 as IClone;

IPoint toPt1=pClone.Clone();

ITransform2D pTrans2D=toPt1 as ITransform2D;

pTrans2D.Move(Dx，Dy);

pTrans2D.Rotate(pPoint，Drotate);

pTrans2D.Scale(pPoint，Dxscale，Dyscale);

toPtsCol.AddPoint(toPt1，ref obj，ref obj);

……

3)根據控制點坐標進行柵格數據重采樣

∥通過IGeoReference接口或IRasterGeometryProc接口進行柵格數據空間轉換

IGeoReference pGR=pRasterLayer as IGeoReference;

pGR.Warp(fromPtsCol，toPtsCol，0);

pGR.Register();

值得注意的是，在柵格數據的空間轉換過程中，由于對數據進行了重采樣，空間轉換后的相鄰兩幅柵格數據邊緣會存在縫隙，無法做到無縫拼接。一種可行的解決方法是在轉換前將柵格數據劃分為具有一定重疊帶的相鄰圖幅，轉換后再將這些圖幅拼接在一起，以消除縫隙對柵格數據的影響。

四、矢柵數據投影變換

1.投影變換方法

在已知矢柵數據空間轉換前后坐標系參數的情況下，可以使用投影變換方法實現基礎地理數據的空間轉換。ArcGIS軟件平臺的坐標系統主要有地理坐標系統(geographic coordinate system)和投影坐標系統(projected coordinate system)。其中，地理坐標系統定義了大地橢球體和大地基準面參數，如橢球體長半徑、短半徑、偏心率等;投影坐標系統則定義了投影參數，如東偽偏移(false easting)、北偽偏移(false northing)、中央經線(central meridian)等。

矢柵數據的投影變換流程如圖3所示，對于矢量數據來說，仍然采用逐要素的變換方式，使用IGeometry接口的Project方法依次更新每個圖元要素實體，實現不同空間參照系之間的相互轉換;對于柵格數據，則可以根據IRaster GeometryPro接口的ProjectFast方法實現投影變換，并采用ISaveAs接口將柵格存儲為Image、Tiff、Grid等數據格式。

2.投影變換功能實現

在進行投影變換之前，需要確定轉換前后的空間參考信息，定義地理坐標系統和投影坐標系統的相關參數，其主要代碼與說明如下:

1)設置空間參考信息

∥定義投影方式為高斯克呂格投影

IProjection pProj=pFact.CreateProjection((int)esriSR_

ProjectionType.esriSRProjection_GaussKruger);

∥定義投影單位、中央經線、東偽偏移、北緯偏移、比例系數、原點緯度等投影參數

IParameter［］pParm=new IParameter［5］;

pParm［0］=pFact.CreateParameter((int)esriSRParameterType.esriSRParameter_FalseEasting);

……

2)實現投影變換

∥若為矢量數據，則依次對每個圖形要素進行重投影處理，并更新圖形實體

IGeometry pGeo=pFeature.Shape as IGeometry;

pGeo.SpatialReference=pSpatialRefFrom;

pGeo.Project((ISpatialReference)pSpatialRefTo);

∥若為柵格數據，則對整個柵格對象賦予投影信息，并進行重采樣

IRasterGeometryProc pRasterGeometryPropc=new Raster-GeometryProcClass();

pRasLayer.SpatialReference=pSpatialRefFrom;

IRaster pRaster=pRasLayer.Raster;

pRasterGeometryPropc.ProjectFast(pSpatialRefTo，0，ref Missing，pRaster);

圖3 投影變換方法流程圖

五、應用實例與誤差分析

筆者基于二次開發組件ArcGIS Engine，采用C#語言設計開發了矢柵基礎地理數據空間轉換系統，界面如圖4所示。對同一區域的矢量和柵格數據采用同一套轉換參數進行空間轉換，將轉換結果疊加顯示，如圖5所示，矢柵數據轉換后具有較好的一致性。

為驗證轉換精度，選取4幅1∶1 000地形圖從地方獨立坐標系轉換到2000國家大地坐標系，將轉換后的圖形數據與實際數據進行對比，抽樣選取25個同名點進行誤差分析，圖6為轉換后隨機采樣點的坐標分量殘差，表1為誤差統計表。由表1可知，X坐標的中誤差為±7 mm，平均值為0 mm;Y坐標的中誤差為±6 mm，平均值為±3 mm。