轉一個關於GDAL座標變換的文章

原帖出處：http://www.gissky.net/Article/647.htm

1. 座標轉化
座標轉化總是一個令人頭疼的問題。不知道你有沒有寫過座標轉換的程式碼，如果寫過就會有體會。那一長串的公式和成群出現的大括弧，真是太可怕了！只有寫過這種東西，你才會對我們生活的這個家園肅然起敬。才會發現原來地球居然如此複雜而顯得如此偉大，人類的數學和地理學家居然是如此地會折騰而顯得如此偉大。敬畏會讓你從此以後不敢亂丟垃圾和隨地吐痰～～呵呵，話說大了。反正我是不敢碰那些符號和公式了。不過既然踏入GIS的門檻，就踏進了座標的深淵，不得不面對它。算了，還是找一個人家已經做出來的輪子去拼自己的車吧，至於輪子有幾根鋼線，什麼材質，怎麼才能搞的比較圓，我就不費神去想了。
有一種叫PROJ的牌子的輪子不知道聽說過沒有。反正開源界現在最好的有關座標轉換的輪子差不多就是它了。對於C/C++人來說，很幸運，因為這個庫可以直接調用。如果是在linux平台下的java人，也很幸運，因為可以用jni來間接調用。如果你是在windows平台下的java人，可能是不幸的。我和我老大折騰了兩天都沒有編譯成功，查了無數的文章，也沒有發現有解決方案（如果你成功了，一定要告訴我哦～～～）。對於業餘的python人來說，前兩天我還以為我們是不幸的，但是今天我覺得上帝還是偏愛python的，讓我找到了這個，還有這個。嗚啦～～～
對於pyproj，我也是“不太熟悉的”:-)，那我們也就看看osr裡面座標轉換是怎麼玩的吧。 先翻譯官方教程，我翻譯不對的地方請回頭看英文原版

1.1. 官方的教程

1.1.1. 簡介
OGRSpatialReference和OGRCoordinateTransformation類提供了用來描繪座標系統（投影和基準面）以及座標系統相互之間轉換的服務。這些服務在OpenGIS座標轉換說明文檔裡有模型，並且有對應的WKT描述。 一些OpenGIS座標系統和服務的背景資料可以在COM的簡單要素(Simple Features for COM)和空間參考系統抽象模型文檔（可以從opengis.org擷取）中找到。GeoTiff投影變換列表也可以輔助地用來理解WKT中的投影公式。EPSG測地學網頁也是一個有用的資源。

1.1.2. 定義一個地理座標系
地理座標系被封裝進了OGRSpatialReference類。有幾種辦法來初始化OGRSpatialReference對象以形成一個合法的座標系統。有兩種主要的座標系統。一種是地理座標(用經緯度表示)。另一種是投影座標(如UTM，用米等單位量度來定位)。
一個地理座標包含基準面(它包含了由長半軸描述的橢球體和反扁率)，本初子午線（一般是格林威治），和一個角度單位（一般是度）。下面的代碼就初始化了一個地理座標系。它提供了圍繞使用者定義名字的地理座標系的所有資訊。

Toggle line numbers Toggle line numbers
   1     OGRSpatialReference oSRS;
   2         oSRS.SetGeogCS( "My geographic coordinate system",
   3                         "WGS_1984",
   4                         "My WGS84 Spheroid",
   5                         SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,
   6                         "Greenwich", 0.0,
   7                         "degree", SRS_UA_DEGREE_CONV );
在這些值中，"My geographic coordinate system", "My WGS84 Spheroid", "Greenwich" 和 "degree" 不是關鍵字。但是被用於顯示給使用者看。但是基準面名"WGS_1984" 經常被用於定義基準面。而且哪些值可以用哪些不行都是有規矩的。不能亂寫。
OGRSpatialReference 已經支援一些標準的座標系統。比如"NAD27", "NAD83", "WGS72" and "WGS84"。要建造它們只要用一個函數SetWellKnownGeogCS().

--------------------------------------------------------------------------------

oSRS.SetWellKnownGeogCS( "WGS84" );

--------------------------------------------------------------------------------

如果EPSG資料庫存在的話，所有EPSG中的地理座標系都可以用GCS編碼來表示。

--------------------------------------------------------------------------------

oSRS.SetWellKnownGeogCS( "EPSG:4326" );

--------------------------------------------------------------------------------

在各種座標系統表達方式中，WKT是個紐帶，通過它，各種表達方式可以互換。一個OGRSpatialReference可以用一個wkt來初始化，或者轉換出wkt表達。

Toggle line numbers Toggle line numbers
   1 char    *pszWKT = NULL;
   2 oSRS.SetWellKnownGeogCS( "WGS84" );
   3 oSRS.exportToWkt( &pszWKT );
   4 printf( "%s/n", pszWKT );
比如這樣一個wkt：

Toggle line numbers Toggle line numbers
   1 GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,
   2 AUTHORITY["EPSG",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG",6326]],
   3 PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["DMSH",0.0174532925199433,
   4 AUTHORITY["EPSG",9108]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG",
   5 4326]]
或者讓它更好看一點吧：

Toggle line numbers Toggle line numbers
   1 GEOGCS["WGS 84",
   2     DATUM["WGS_1984",
   3         SPHEROID["WGS 84",6378137,298.257223563,
   4             AUTHORITY["EPSG",7030]],
   5         TOWGS84[0,0,0,0,0,0,0],
   6         AUTHORITY["EPSG",6326]],
   7     PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
   8     UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],
   9     AXIS["Lat",NORTH],
  10     AXIS["Long",EAST],
  11     AUTHORITY["EPSG",4326]]
OGRSpatialReference::importFromWkt()方法可以用來把wkt座標系統設定到OGRSpatialReference中。

 

1.1.3. 定義一個投影系統
一個投影座標系統（如UTM等）基於一個地理座標系統定義以便於線上性單位和角度位置之間轉換。接下來的代碼定義了一個在WGS84基準面下的地理座標系統為基礎的UTM17帶投影座標系統。

Toggle line numbers Toggle line numbers
   1         OGRSpatialReference     oSRS;
   2         oSRS.SetProjCS( "UTM 17 (WGS84) in northern hemisphere." );
   3         oSRS.SetWellKnownGeogCS( "WGS84" );
   4         oSRS.SetUTM( 17, TRUE );
調用SetProjCS()設定一個使用者定義名字的投影座標系統並確定系統被投影過。SetWellKnownGeogCS()分配一個地理座標系統，SetUTM()設定投影轉換參數細節。上面的步驟直到這時才能建立一個合法的定義。但是這個對象需要時將會重新自動編排內在表達，以保持合法性。
請尤其注意定義對象時的步驟順序！！
上面的定義會定義一個像下面一樣的WKT。注意UTM117被擴充成橫軸莫卡托定義UTM帶的參數。

Toggle line numbers Toggle line numbers
   1 PROJCS["UTM 17 (WGS84) in northern hemisphere.",
   2     GEOGCS["WGS 84",
   3         DATUM["WGS_1984",
   4             SPHEROID["WGS 84",6378137,298.257223563,
   5                 AUTHORITY["EPSG",7030]],
   6             TOWGS84[0,0,0,0,0,0,0],
   7             AUTHORITY["EPSG",6326]],
   8         PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
   9         UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],
  10         AXIS["Lat",NORTH],
  11         AXIS["Long",EAST],
  12         AUTHORITY["EPSG",4326]],
  13     PROJECTION["Transverse_Mercator"],
  14     PARAMETER["latitude_of_origin",0],
  15     PARAMETER["central_meridian",-81],
  16     PARAMETER["scale_factor",0.9996],
  17     PARAMETER["false_easting",500000],
  18     PARAMETER["false_northing",0]]
不同的投影會有不同的函數來設定參數(SetTM() (Transverse Mercator), SetLCC() (Lambert Conformal Conic)和 SetMercator())

1.1.4. 查詢座標系統
一旦OGRSpatialReference建立起來了，就可以查詢多種資訊。用OGRSpatialReference::IsProjected() 和 OGRSpatialReference::IsGeographic()可以看他是否是投影座標還是地理座標。OGRSpatialReference::GetSemiMajor(), OGRSpatialReference::GetSemiMinor() 和OGRSpatialReference::GetInvFlattening()可以擷取橢球體資訊。OGRSpatialReference::GetAttrValue() 可以用來擷取PROJCS, GEOGCS, DATUM, SPHEROID, 和PROJECTION 的名字表達字串。OGRSpatialReference::GetProjParm()可以用來擷取投影參數。OGRSpatialReference::GetLinearUnits() 可以用來擷取線性單位，並轉換成米。
下面的代碼示範了GetAttrValue()擷取投影的用法，以及GetProjParm()GetAttrValue() 的用法。已經定義的投影參數(如SRS_PP_CENTRAL_MERIDIAN)應該在GetProjParm()擷取投影參數的時候使用。ogrspatialreference.cpp中設定不同投影的方法的代碼可以用來尋找哪個投影用哪個參數。

Toggle line numbers Toggle line numbers
   1 const char *pszProjection = poSRS->GetAttrValue("PROJECTION");
   2     if( pszProjection == NULL )
   3     {
   4         if( poSRS->IsGeographic() )
   5             sprintf( szProj4+strlen(szProj4), "+proj=longlat " );
   6         else
   7             sprintf( szProj4+strlen(szProj4), "unknown " );
   8     }
   9     else if( EQUAL(pszProjection,SRS_PT_CYLINDRICAL_EQUAL_AREA) )
  10     {
  11         sprintf( szProj4+strlen(szProj4),
  12            "+proj=cea +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f ",
  13                  poSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
  14                  poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
  15                  poSRS->GetProjParm(SRS_PP_FALSE_EASTING,0.0),
  16                  poSRS->GetProjParm(SRS_PP_FALSE_NORTHING,0.0) );
  17     }
  18     ...

1.1.5. 座標轉換
OGRCoordinateTransformation 類用來在不同座標系統間轉換座標。新的轉換對象由OGRCoordinateTransformation()建立。建立後OGRCoordinateTransformation::Transform() 方法就可以用來轉換不同座標系的點了。

Toggle line numbers Toggle line numbers
   1         OGRSpatialReference oSourceSRS, oTargetSRS;
   2         OGRCoordinateTransformation *poCT;
   3         double                  x, y;
   4         oSourceSRS.importFromEPSG( atoi(papszArgv[i+1]) );
   5         oTargetSRS.importFromEPSG( atoi(papszArgv[i+2]) );
   6         poCT = OGRCreateCoordinateTransformation( &oSourceSRS,
   7                                                   &oTargetSRS );
   8         x = atof( papszArgv[i+3] );
   9         y = atof( papszArgv[i+4] );
  10         if( poCT == NULL || !poCT->Transform( 1, &x, &y ) )
  11             printf( "Transformation failed./n" );
  12         else
  13             printf( "(%f,%f) -> (%f,%f)/n",
  14                     atof( papszArgv[i+3] ),
  15                     atof( papszArgv[i+4] ),
  16                     x, y );
有兩個原因導致在轉換時會出錯。第一，OGRCreateCoordinateTransformation() 可能失敗。一般是因為程式認為兩個指定座標系統不能建立轉換關係，這可能是由於使用的投影Proj內部暫時不支援，這樣兩個不一樣的基準面就沒有已知關係。或者一個投影定義得不適當。如果OGRCreateCoordinateTransformation() 失敗了將返回空。
OGRCoordinateTransformation::Transform() 方法本身也可能失敗。這可能是因為一個有上面問題積累後形成的問題。或者是一個“因為一個或者多個點因數值上未定義而省略操作”起因的後果。Transform()函數如果成功，返回TRUE，如果某個點轉換失敗，剩下的點就會以一個錯誤的狀態保持。
座標轉換服務其實可以處理3D點的。並且服務會根據橢球體和基準面上的高程差異調節高程。將來，也有可能會有在不同向量基準面上的點間轉換的應用。如果沒有Z值，則假設點都在大地水準面上。
下面的例子顯示了如何方便地建立一個地理座標系統到一個基於同一個地理座標系統的投影座標系統。並在兩個座標系統之間轉換。

Toggle line numbers Toggle line numbers
   1     OGRSpatialReference    oUTM, *poLatLong;
   2     OGRCoordinateTransformation *poTransform;
   3     oUTM.SetProjCS("UTM 17 / WGS84");
   4     oUTM.SetWellKnownGeogCS( "WGS84" );
   5     oUTM.SetUTM( 17 );
   6     poLatLong = oUTM.CloneGeogCS();
   7     poTransform = OGRCreateCoordinateTransformation( &oUTM, poLatLong );
   8     if( poTransform == NULL )
   9     {
  10         ...
  11     }
  12     ...
  13     if( !poTransform->Transform( nPoints, x, y, z ) )
  14     ...

1.1.6. API綁定
C的介面在ogr_srs_api.h中定義，python綁定在osr.py模組中。一些C++的方法在C和python綁定中沒有定義。
Python Bindings(C綁定略過)

Toggle line numbers Toggle line numbers
   1 class osr.SpatialReference
   2     def __init__(self,obj=None):
   3     def ImportFromWkt( self, wkt ):
   4     def ExportToWkt(self):
   5     def ImportFromEPSG(self,code):
   6     def IsGeographic(self):
   7     def IsProjected(self):
   8     def GetAttrValue(self, name, child = 0):
   9     def SetAttrValue(self, name, value):
  10     def SetWellKnownGeogCS(self, name):
  11     def SetProjCS(self, name = "unnamed" ):
  12     def IsSameGeogCS(self, other):
  13     def IsSame(self, other):
  14     def SetLinearUnits(self, units_name, to_meters ):
  15     def SetUTM(self, zone, is_north = 1):
  16 class CoordinateTransformation:
  17     def __init__(self,source,target):
  18     def TransformPoint(self, x, y, z = 0):
  19     def TransformPoints(self, points):

1.1.7. 介面實現
OGRCoordinateTransformation服務的實現是建立在PROJ4庫之上，該庫是最先由USGS的Gerald Evenden 編寫的。

 

1.2. 我的實驗

1.2.1. 建立和比較

Toggle line numbers Toggle line numbers
   1 >>> import gdal
   2 >>> import osr
   3 >>> dataset = gdal.Open("e:/gisdata/gtif/spot.tif")
從資料集中擷取空間參考並且建立一個SpatialReference對象

Toggle line numbers Toggle line numbers
   1 >>> sr = dataset.GetProjectionRef()
   2 >>> sr
   3 'PROJCS["unnamed",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Cla
   4 rke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG"
   5 ,"6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPS
   6 G","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
   7 PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["f
   8 alse_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EP
   9 SG","9001"]],AUTHORITY["EPSG","26713"]]'
  10 >>> osrobj = osr.SpatialReference()
  11 >>> osrobj.ImportFromWkt(sr)
  12 0
  13 >>> osrobj.ExportToWkt()
  14 'PROJCS["unnamed",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Cla
  15 rke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG"
  16 ,"6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPS
  17 G","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
  18 PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["f
  19 alse_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EP
  20 SG","9001"]],AUTHORITY["EPSG","26713"]]'
  21 >>> osrobj.IsGeographic()
  22 0
  23 >>> osrobj.IsProjected()
  24 1
ERROR: EOF in multi-line statement

再擷取一個進行比較

Toggle line numbers Toggle line numbers
   1 >>> dataset2 = gdal.Open("e:/gisdata/gtif/uparea.tif")
   2 >>> sr2 = dataset2.GetProjectionRef()
   3 >>> sr2
   4 'PROJCS["unnamed",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Cla
   5 rke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG"
   6 ,"6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPS
   7 G","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
   8 PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["f
   9 alse_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EP
  10 SG","9001"]],AUTHORITY["EPSG","26713"]]'
  11 >>> osrobj2 = osr.SpatialReference()
  12 >>> osrobj2.ImportFromWkt(sr2)
  13 0
  14 >>> osrobj2.IsSame(osrobj)
  15 1
ERROR: EOF in multi-line statement

建立一個地理座標系，然後和已有的座標系統比較

Toggle line numbers Toggle line numbers
   1 >>> osrobj3 = osr.SpatialReference()
   2 >>> osrobj3.SetWellKnownGeogCS("WGS84")
   3 0
   4 >>> osrobj3.IsSame(osrobj2)
   5 0
   6 >>> osrobj3.IsSame(osrobj)
   7 0
   8 >>> osrobj3.ExportToWkt()
   9 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHOR
  10 ITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Gre
  11 enwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["
  12 EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]]'
  13 >>> osrobj3.IsGeographic()
  14 1

1.2.2. 地理座標系和投影座標系之間的座標轉換
spot的映像投影座標範圍是這樣的

Toggle line numbers Toggle line numbers
   1 In projected or local coordinates
   2 Left: 590000.000000
   3 Right: 609000.000000
   4 Top: 4928000.000000
   5 Bottom: 4914000.000000
用CoordinateTransformation轉換結果:

Toggle line numbers Toggle line numbers
   1 >>> ct = osr.CoordinateTransformation(osrobj,osrobj3)
   2 >>> ct.TransformPoint(590000,4928000)
   3 (-103.86789483706976, 44.501658895816213, 2.1979212760925293e-007)
   4 >>> ct.TransformPoint(609000,4928000)
   5 (-103.6289570318164, 44.499040134967487, 2.2072345018386841e-007)
   6 >>> ct.TransformPoint(609000,4914000)
   7 (-103.63190060724861, 44.373035580280714, 2.2072345018386841e-007)
   8 >>> ct.TransformPoint(590000,4914000)
   9 (-103.87032571673741, 44.375642923223602, 2.1979212760925293e-007)
  10 >>>
ArcGIS裡面的地理座標範圍結果是這樣的：

Toggle line numbers Toggle line numbers
   1 In decimal degrees
   2 West: -103.870326
   3 East: -103.628957
   4 North: 44.501659
   5 South: 44.373036
用proj計算的和ArcGIS有些差距(分別有一個很准，一個差了0.01度)。但是我覺得這種誤差應該是允許的。