GIS空間參考及坐標轉換

 

 

空間參考（Spatial Reference）是 GIS 數據的骨骼框架，能夠將我們的數據定位到相應的位置，為地圖中的每一點提供准確的坐標。 在同一個地圖上顯示的地圖數據的空間參考必須是一致的，如果兩個圖層的空間參考不一致，往往會導致兩幅地圖無法正確拼合，因此開發一個 GIS 系統時，為數據選擇正確的空間參考非常重要。

1 相關知識 

1.1大地水准面 

大地水准面是由靜止海水面並向大陸延伸所形成的不規則的封閉曲面。

1.2地球橢球體 

由定義可以知大地水准面的形狀也是不規則的，仍不能用簡單的數學公式表示，為了測量成果的 計算和制圖的需要，人們選用一個同大地水准面相近的可以用數學方法來表達的橢球體來代替， 簡稱地球橢球體，它是一個規則的曲面，是測量和制圖的基礎，因地球橢球體是人們選定的跟大 地水准面很接近的規則的曲面，所以地球橢球體就可以有多個，地球橢球體是用長半軸、短半軸 和扁率來表示的。

橢球名稱	長半軸（m）	短半軸（m）	扁率的倒數，1/f
克拉克 (Clarke 1866)	6378 206.4	6356 583.8	294.978 8 2
貝塞爾 （Bessel 1841）	6 377 397.155	6 356 078.965	299.152 843 4
International 1924	6 378 388	6 356 911.9	296.999 362 1
克拉索夫斯基 （1940）	6 378 245	6 356 863	298.299 738 1
GRS 1980	6 378 137	6 356 752.3141	298.257 222 101
WGS 1984	6 378 137	6 356 752.3142	298.257 223 563

1.3基准面 

 

基准面是在特定區域內與地球表面極為吻合的橢球體。橢球體表面上的點與地球表面上的特定位 置相匹配,也就是對橢球體進行定位,該點也被稱作基准面的原點。原點的坐標是固定的，所有其他點由其計算獲得

基准面的坐標系原點往往距地心有一定偏移(有的也在地心,如 WGS1984), 如西安 80 的基准面和北京 54的基准面。因為橢球體通過定位以便能更好的擬合不同的地區,所以同一個橢球體可以擬合好幾個基准面.因為原點不同,所以不同的基准面上,同一個點的坐標是不相同的,這點我們應該清楚.下面以華盛頓州貝靈厄姆市為例來說明。使用 NAD27、NAD83 和 WGS84 以十進制為單位比較貝靈厄姆的坐標。顯而易見，NAD83 和 WGS84 表示的坐標幾乎相同，但 NAD27 表示的坐標則大不相同，這是因為所使用的基准面和旋轉橢球體對地球基本形狀的表示方式不同。

1.4 地圖投影 

簡單的說地圖投影就是把地球表面的任意點，利用一定數學法則，轉換到地圖平面上的理論和方法.2.

 

2. 兩種坐標系 

2.1 地理坐標系 

地理坐標系也可稱為真實世界的坐標系，是用於確定地物在地球上位置的坐標系，它用經緯度來表示地物的位置，經度和緯度是從地心到地球表面上某點的測量角，通常以度或百分度為單位來測量該角度。

2.2 投影坐標系 

投影坐標系是基於地理坐標系的，它使用基於 X,Y 值的坐標系統來描述地球上某個點所處的位置，可以這樣認為投影坐標系=地理坐標系(如：北京 54、西安 80、WGS84)+投影方法(如：高斯－克呂格、Lambert 影、Mercator 投影)+線性單位。

2.3 ArcGIS Engine 對空間參考的支持 

ArcGIS Engine 提供了一系列對象供開發者管理 GIS 系統的坐標系統。對大部分開發者而言了解ProjectedCoordinateSystem, GeographicCoordinateSystem, SpatialReference Environment 這三個組件類是非常有必要的，對於高級開發者而言，可能需要自定義坐標系統可以使用這些對象Projection,Datum,AngularUnit,Spheriod,PrimeMeridian 和 GeoTransformation 等。在 ArcGIS 中除了我們上面介紹的兩種坐標系，還有一種稱之為 Unknown 的坐標系，這種坐標系是當我們的數據沒有坐標（jpg等文件）或者坐標文件丟失的時候 ArcMap 不能識別數據的投影信息而賦予的，在 ArcGIS Engine 中下面三個類分別對應了三個坐標系： 

 

利用 ArcGIS Engine 創 建 一 個 坐 標 系 或 者 基准面用的是SpatialReferenceEnvironmentClass 類，該類實現了 ISpatialReferenceFactory 接口，該接口 定義了創建坐標系，基准面等方法和屬性。

在利用 ISpatialReferenceFactory 創建坐標系的時候往往需要一個 int 類型的參數，這個int 其實就是這些坐標系的代號，如我們熟悉的 4326 就是 WGS1984。

2.4 同一基准面的坐標轉換 

對於同一基准面，我們可以肯定一點就是同一位置經緯度坐標是一樣的，而不同的就是計算成平面坐標的時候可能有所不同，因為算法不一樣，在這里我只是將經緯度坐標轉成平面的坐標。

// 同一橢球下經緯度坐標、平面的坐標互轉
private IPoint GetpProjectPoint(IPoint pPoint, bool pBool)         
 {             
     ISpatialReferenceFactory pSpatialReferenceEnvironemnt = new SpatialReferenceEnvironment();             
     ISpatialReference pFromSpatialReference = pSpatialReferenceEnvironemnt.CreateGeographicCoordinateSystem((int)esriSRGeoCS3Type.esriSRGeoCS_Xian1980);//西安80 
     ISpatialReference pToSpatialReference = pSpatialReferenceEnvironemnt.CreateProjectedCoordinateSystem((int)esriSRProjCS4Type.esriSRProjCS_Xian1980_3_Degree_GK_Zone_34);//西安80             
     if (pBool == true)//球面轉平面            
        {                   
            Geometry pGeo = (IGeometry)pPoint;                 
            pGeo.SpatialReference = pFromSpatialReference;               
            pGeo.Project(pToSpatialReference);                 
            return pPoint;             
        }             
        else //平面轉球面             
        {  
        IGeometry pGeo = (IGeometry)pPoint;                 
        pGeo.SpatialReference = pToSpatialReference;                
        pGeo.Project(pFromSpatialReference);                 
        return pPoint;             
        }         
}

2.5 不同基准面的坐標轉換 

2.5.1 三參數與七參數

通過前面的介紹，我們知道地球上同一位置的坐標在不同的基准面上是不一樣的，而基准面是構成坐標系的一個部分，因為基准面在定位的時候牽扯到了相對地心的平移或旋轉等，所以對於這樣的轉換我們無法直接進行，需要一個轉換參數，而這些參數也是基於不同的模型的，常用的有三參數和 7 參數,三參數是比較簡單的也是比較容易理解的，三參數是在兩個基准面之間進行了 X, Y, Z 軸的平移，通過下面的圖我們很清楚的看到三參數之間兩個基准面的關系：

如果知道了這三個平移的參數（dx, dy, dz），外加1個基准面上的點，那么：

而 7 參數的模型比較復雜，這種復雜的同時讓精度大為提高，7參數不僅僅考慮了兩個基准面之間的平移，還考慮了旋轉外加一個比例因子(橢球體的大小可能不一樣).

空間直角坐標系(XYZ) 轉換(To) 空間直角坐標系(XYZ)：（布爾沙模型，此步重要，將一個橢球基准轉換到另一個橢球基准）

其中△X，△Y，△Z是坐標平移量，R（ω）是旋轉矩陣，（1+m）是比例因子。

7參數轉換整體流程（不同橢球之間）：

平面直角坐標系-->大地坐標系-->空間直角坐標系      Convert to      空間直角坐標系-->大地坐標系-->平面直角坐標系

舉個栗子，比如從BJ1954平面直角坐標系 轉到WGS84平面直角坐標系那么需要5步：

①BJ1954平面直角坐標系 至 BJ1954大地坐標系

②BJ1954大地坐標系 至BJ1954空間直角坐標系

③BJ1954空間直角坐標系 至WGS84空間直角坐標系

④WGS84空間直角坐標系 至WGS84大地坐標系

⑤WGS84大地坐標系 至WGS84平面直角坐標系

如果從BJ1954空間直角坐標系轉到WGS84平面直角坐標系只需一步：

①BJ1954空間直角坐標系 至 WGS84空間直角坐標系

2.5.2 最小二乘法求解七、四參數

七參數求解大致流程為：

a)

b) 如上5-7公式符合最小二乘法A*X=B的形式，依照最小二乘法求解：

四參數轉換整體流程：（同一橢球下不同坐標系,注意與七參數的區別）

大地坐標系-->平面直角坐標系-->平面直角坐標系-->大地坐標系

 

 

2.5.3 DEMO坐標轉換的整體流程

 

 

2.5.4 七參數的ArcGIS編程實現

對於不同基准面之間的轉換，ArcGIS Engine 提供了一個用來控制轉換參數的接口 IGeoTransformation，該接口被被較多的其它類實現着。

每一個接口對應了一種轉換方法，比如 GeocentricTranslationClass類就實現了三參數，而CoordinateFrameTransformationClass 類實現了7 參數，要實現 3 參數或者 7 參數需要 IGeometry2 或更新接口的ProjectEx 方法，下面我們用代碼實現一個不同基准面之間的坐標轉換。 

 

// 不同基准面之間的坐標轉換 WGS1984 GaussKruger -- > esriSRProjCS_Xian1980_GK_Zone_19
public void ProjectExExample()        
{                          
    ISpatialReferenceFactory pSpatialReferenceFactory = new SpatialReferenceEnvironmentClass();  // ISpatialReference pFromCustom = pSpatialReferenceFactory.CreateESRISpatialReferenceFromPRJFile(@"E:\arcgis\Engine\z idingyi.prj");  
    
    IPoint pFromPoint = new PointClass();  
    pFromPoint.X = 518950.788;  
    pFromPoint.Y = 4335923.97;  
    IZAware pZAware = pFromPoint as IZAware;  
    pZAware.ZAware = true;  
    pFromPoint.Z = 958.4791;  // 自定義投影WGS84下的北京6度 19帶。             
    ((IGeometry)pFromPoint).SpatialReference = CreateCustomProjectedCoordinateSystem();  // ((IGeometry)pFromPoint).SpatialReference = pFromCustom; 
         
    //目標投影             
    IProjectedCoordinateSystem projectedCoordinateSystem = pSpatialReferenceFactory.CreateProjectedCoordinateSystem((int)esriSRProjCS4Type.esriSRProjCS_Xian1980_GK_Zone_19);             //因為目標基准面和原始基准面不在同一個上，所以牽扯到參數轉換，我用7參數 轉換             ICoordinateFrameTransformation   pCoordinateFrameTransformation = new CoordinateFrameTransformationClass();  
    pCoordinateFrameTransformation.PutParameters(-112.117, 4.530, 21.89, -0.00058702, -0.00476421, 0.00009358, 0.99998006411);   
    pCoordinateFrameTransformation.PutSpatialReferences(CreateCustomProjectedCoordinateSystem(), projectedCoordinateSystem as ISpatialReference);                       

    //投影轉換             
    IGeometry2 pGeometry = pFromPoint as IGeometry2;  
    pGeometry.ProjectEx(projectedCoordinateSystem as ISpatialReference, esriTransformDirection.esriTransformForward, pCoordinateFrameTransformation, false, 0, 0);  
}

private IProjectedCoordinateSystem CreateCustomProjectedCoordinateSystem()
{
    ISpatialReferenceFactory2 pSpatialReferenceFactory = new SpatialReferenceEnvironmentClass();
    IProjectionGEN pProjection = pSpatialReferenceFactory.CreateProjection((int)esriSRProjectionType.esriSRProjection_GaussKruger) as IProjectionGEN; 
    IGeographicCoordinateSystem pGeographicCoordinateSystem = pSpatialReferenceFactory.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);
    ILinearUnit pUnit = pSpatialReferenceFactory.CreateUnit((int)esriSRUnitType.esriSRUnit_Meter) as ILinearUnit;
    IParameter[] pParameters = pProjection.GetDefaultParameters();

    IProjectedCoordinateSystemEdit pProjectedCoordinateSystemEdit = new ProjectedCoordinateSystemClass();
    object pName = "WGS-BeiJing1954";
    object pAlias = "WGS-BeiJing1954";
    object pAbbreviation = "WGS-BeiJing1954";
    object pRemarks = "WGS-BeiJing1954";
    object pUsage = "Calculate Meter From lat and lon";
    object pGeographicCoordinateSystemObject = pGeographicCoordinateSystem as object;
    object pUnitObject = pUnit as object;
    object pProjectionObject = pProjection as object;
    object pParametersObject = pParameters as object;
    pProjectedCoordinateSystemEdit.Define(ref pName, ref pAlias, ref pAbbreviation, ref  pRemarks, ref pUsage, ref pGeographicCoordinateSystemObject, 
        ref pUnitObject, ref pProjectionObject, ref pParametersObject);
    IProjectedCoordinateSystem5 pProjectedCoordinateSystem = pProjectedCoordinateSystemEdit as IProjectedCoordinateSystem5;
    pProjectedCoordinateSystem.FalseEasting = 500000;
    pProjectedCoordinateSystem.LatitudeOfOrigin = 0;
    pProjectedCoordinateSystem.set_CentralMeridian(true, 111);
    pProjectedCoordinateSystem.ScaleFactor = 1;
    pProjectedCoordinateSystem.FalseNorthing = 0;
    
    return pProjectedCoordinateSystem;
}

 

 

參考文章

坐標轉換流程與公式 七參數 四參數，2017.8.15.

ArcGIS E0ngine開發之旅10--空間參考及坐標轉換，2016.8.4

在ArcGIS Desktop中進行三參數或七參數精確投影轉換

墨卡托投影、高斯-克呂格投影、UTM投影及我國分帶方法

地圖投影與ArcGIS坐標系轉換

地理坐標系、大地坐標系與地圖投影與重投影詳解

安衛. 一種平面四參數法坐標轉換方法的實現, 2012.