GDAL投影轉換

轉自 http://blog.csdn.net/jingss_3/article/details/8240328

一、簡介

本文參考英文地址：http://www.gdal.org/ogr/osr_tutorial.html。

OGRSpatialReference類和OGRCoordinateTransformation類主要用來提供定義坐標系統（投影和水准面）和轉換坐標。這兩個類都基於OpenGIS的坐標轉換說明，並且使用Well Known Text格式來進行表述坐標系統。

一些關於OpenGIS坐標系統的資料，以及空間參考坐標抽象模型文件可以在OGC（Open Geospatial Consortium）的網站上找到。GeoTIFF投影轉換列表（GeoTIFF Projections Transform List）網頁可以更好的幫助你理解WKT的規則，同時EPSG的網站也是很有用的資料。

二、定義地理坐標系統

坐標系統使用OGRSpatialReference類來進行封裝。這里提供了數種初始化OGRSpatialReference類的方式。這里有兩類主要的坐標系統，第一種是地理坐標系統（位置信息使用經緯度來表示的），第二種是投影坐標系統（比如UTM-通用橫軸墨卡托投影，位置信息使用米或者英尺來表示）。

一個地理坐標系統包含的信息有一個大地基准（里面含有一個使用長半軸和扁率的倒數來表示的托球體），一個中央經線（通常是本初子午線，也就是0度經線Greenwich）, 此外還有一個角度的度量單位，使用度而不是弧度。如果含有這些信息，就可以構造一個有效的地理坐標系統。

OGRSpatialReference oSRS;
oSRS.SetGeogCS( "Mygeographic coordinate system",
               "WGS_1984",
               "My WGS84 Spheroid",
               SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,
               "Greenwich", 0.0,
                "degree", SRS_UA_DEGREE_CONV );

在上面的代碼中，名稱為“My geographic coordinate system”，“My WGS84 Spheroid”，“Greenwich”和“degree”的並不是關鍵詞，這些主要是用來給用戶進行說明的。然而，“WGS_1984”是一個定義大地基准的關鍵詞，注意：這里的大地基准必須是一個有效的大地基准！（這句話的意思，前面的那些字符串就是隨便指定的，用來顯示的，后面的WGS_1984這個位置的字符串，必須是一個有效的，不能隨便命名，具體后面會說到）。

OGRSpatialReference可以使用一些常用的字符串來進行建立一個常用的坐標系統，這些字符串包括“NAD27”、“NAD83”，“WGS72”和“WGS84”等，使用的函數是SetWellKnownGeogCS()，使用方式見下面。

oSRS.SetWellKnownGeogCS( "WGS84" );

而且，還可以使用EPSG數據庫定義的GCS代碼來定義一個地理坐標系統，如：

oSRS.SetWellKnownGeogCS( "EPSG:4326" );

為了方便的和其他庫進行關聯，這里可以轉換為OpenGIS的WKT格式。同樣OGRSpatialReference可以使用一個WKT來進行初始化，或者將里面的信息導出為WKT格式。

char    *pszWKT =NULL;
SRS.SetWellKnownGeogCS( "WGS84" );
oSRS.exportToWkt( &pszWKT );
printf( "%s\n", pszWKT );

上面的語句會輸出下面的內容：

GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG",6326]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG",4326]]

將上面的字符串格式調整成更好看的樣子：

GEOGCS["WGS 84",
   DATUM["WGS_1984",
       SPHEROID["WGS 84",6378137,298.257223563,
           AUTHORITY["EPSG",7030]],
       TOWGS84[0,0,0,0,0,0,0],
       AUTHORITY["EPSG",6326]],
    PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
   UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],
   AXIS["Lat",NORTH],
   AXIS["Long",EAST],
   AUTHORITY["EPSG",4326]]

函數OGRSpatialReference::importFromWkt()可以從一個WKT定義的坐標系統來構造一個OGRSpatialReference類對象。

三、定義投影坐標系統

一個投影坐標系統（比如UTM，蘭伯特等角圓錐投影等）需要建立在一個地理坐標系統之上，在投影坐標系統中，坐標點使用米或者英尺等長度單位來表示，同時也可以用經緯度的角度坐標來表示。下面將定義一個UTM的第17帶的投影坐標系統，基於WGS84的大地基准橢球體。

OGRSpatialReference    oSRS;
oSRS.SetProjCS( "UTM 17(WGS84) in northern hemisphere." );
oSRS.SetWellKnownGeogCS("WGS84" );
oSRS.SetUTM( 17, TRUE );

首先調用SetProjCS()函數設置投影坐標系統的名稱，然后使用函數SetWellKnownGeogCS()指定地理坐標系統，最后調用函數SetUTM()設置投影轉換參數信息。完成這些工作之后就定義了一個有效的投影坐標系統。這里必須要注意定義OGRSpatialReference的順序！

上面定義的投影使用WKT表示的形式如下，注意UTM17會使用橫軸墨卡托的分帶投影參數來表示。

PROJCS["UTM 17 (WGS84) in northernhemisphere.",
   GEOGCS["WGS 84",
       DATUM["WGS_1984",
           SPHEROID["WGS 84",6378137,298.257223563,
               AUTHORITY["EPSG",7030]],
           TOWGS84[0,0,0,0,0,0,0],
           AUTHORITY["EPSG",6326]],
        PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
       UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],
       AXIS["Lat",NORTH],
       AXIS["Long",EAST],
       AUTHORITY["EPSG",4326]],
   PROJECTION["Transverse_Mercator"],
   PARAMETER["latitude_of_origin",0],
   PARAMETER["central_meridian",-81],
   PARAMETER["scale_factor",0.9996],
   PARAMETER["false_easting",500000],
   PARAMETER["false_northing",0]]

這里提供了很多設置投影坐標的方法，包括SetTM()(橫軸墨卡托投影), SetLCC()(蘭伯特等角圓錐投影)和SetMercator()。

四、獲取坐標系統信息

一旦一個OGRSpatialReference對象進行創建，那么就可以獲取里面的各種信息，比如可以使用OGRSpatialReference::IsProjected()和OGRSpatialReference::IsGeographic()方法來判斷坐標系統是地理坐標系統還是投影坐標系統。使用函數OGRSpatialReference::GetSemiMajor()、OGRSpatialReference::GetSemiMinor()和OGRSpatialReference::GetInvFlattening()方法可以用來獲取橢球體信息，分別是橢球體的長半軸，短半軸以及扁率的倒數。使用OGRSpatialReference::GetAttrValue()方法可以用來獲取PROJCS、GEOGCS、DATUM、SPHEROID和PROJECTION的名稱字符串。使用OGRSpatialReference::GetProjParm()方法可以獲取投影參數信息。使用OGRSpatialReference::GetLinearUnits()方法可以獲取長度單位類型，並將其轉換為米。

下面的代碼是一個獲取坐標系統信息的示例，摘自ogr_srs_proj4.cpp文件中。使用GetAttrValue()獲取投影名稱，使用GetProjParm()獲取投影參數信息。GetAttrValue()函數的第一個值節點就是WKT字符串的描述信息。投影參數的宏定義（比如SRS_PP_CENTRAL_MERIDIAN）和函數GetProjParm()一起使用，可以用來獲取投影的參數。更多的使用方式可以參考文件ogrspatialreference.cpp中的相關代碼。

const char *pszProjection = poSRS->GetAttrValue("PROJECTION");
if( pszProjection == NULL )
{
    if( poSRS->IsGeographic() )
        sprintf(szProj4+strlen(szProj4), "+proj=longlat " );
    else
        sprintf(szProj4+strlen(szProj4), "unknown " );
}
else if(EQUAL(pszProjection,SRS_PT_CYLINDRICAL_EQUAL_AREA) )
{
    sprintf(szProj4+strlen(szProj4),
      "+proj=cea +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f ",
           poSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
           poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
           poSRS->GetProjParm(SRS_PP_FALSE_EASTING,0.0),
           poSRS->GetProjParm(SRS_PP_FALSE_NORTHING,0.0));
}
...
const char *pszProjection = poSRS->GetAttrValue("PROJECTION");
if( pszProjection == NULL )
{
    if( poSRS->IsGeographic() )
        sprintf(szProj4+strlen(szProj4), "+proj=longlat " );
    else
        sprintf(szProj4+strlen(szProj4), "unknown " );
}
else if(EQUAL(pszProjection,SRS_PT_CYLINDRICAL_EQUAL_AREA) )
{
    sprintf(szProj4+strlen(szProj4),
      "+proj=cea +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f ",
           poSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
           poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
           poSRS->GetProjParm(SRS_PP_FALSE_EASTING,0.0),
           poSRS->GetProjParm(SRS_PP_FALSE_NORTHING,0.0));
}
...

五、坐標轉換

OGRCoordinateTransformation類可以用來在不同的坐標系統中進行坐標轉換。可以使用函數OGRCreateCoordinateTransformation()創建一個新的坐標轉換對象，然后使用OGRCoordinateTransformation::Transform()方法來進行坐標轉換。

OGRSpatialReference oSourceSRS,oTargetSRS;
OGRCoordinateTransformation *poCT;
double                 x, y;

oSourceSRS.importFromEPSG(atoi(papszArgv[i+1]) );
oTargetSRS.importFromEPSG(atoi(papszArgv[i+2]) );

poCT = OGRCreateCoordinateTransformation(&oSourceSRS,
                                         &oTargetSRS );
x = atof( papszArgv[i+3] );
y = atof( papszArgv[i+4] );

if( poCT == NULL || !poCT->Transform( 1, &x,&y ) )
    printf("Transformation failed.\n" );
else
    printf("(%f,%f) -> (%f,%f)\n",
            atof(papszArgv[i+3] ),
            atof(papszArgv[i+4] ),
            x, y );

一對點轉換失敗的原因有，首先，OGRCreateCoordinateTransformation()函數執行失敗，通常的原因是不知道指定的兩個投影之間的轉換關系。這個可能是試用了PROJ.4庫不支持的投影，不同的橢球體之間的轉換關系沒有定義，或者是其中的一個坐標系統沒有定義完全。如果函數OGRCreateCoordinateTransformation() 執行失敗，那么其返回值將是NULL。

OGRCoordinateTransformation::Transform() 方法本身頁可能執行失敗。這個可能的原因也有上面的問題，或者是轉換的點里面有一個以上沒有定義的數字。函數Transform()執行成功是返回TRUE，如果有一個點轉換失敗都會返回FALSE。

坐標轉換也可以支持三維點的坐標轉換，會自動根據不同的橢球地的基准面調整高程值。這個可以用來在不同的垂直基准面進行坐標轉換。如果沒有Z值，系統會認為所有的點都是在水准面上。

下面的例子演示了如果從一個投影坐標系統中的點轉換為該投影中的地理坐標系統中的點，將米表示的坐標轉換為經緯度表示的坐標。

OGRSpatialReference   oUTM, *poLatLong;
OGRCoordinateTransformation *poTransform;
oUTM.SetProjCS("UTM 17 /WGS84");
oUTM.SetWellKnownGeogCS("WGS84" );
oUTM.SetUTM( 17 );
poLatLong = oUTM.CloneGeogCS();

poTransform = OGRCreateCoordinateTransformation( &oUTM,poLatLong );
if( poTransform == NULL )
{
    ...
}

...
if( !poTransform->Transform( nPoints, x,y, z ) )
...

六、其他語言接口

OGR的空間參考提供了一個C語言的接口，定義在ogr_srs_api.h文件中，Python的接口定義在osr.py文件中。所有的接口名稱和C++的接口都很相似，但是C和Python中有些方法沒有進行提供。

1、C語言接口

typedef void *OGRSpatialReferenceH;
typedef void *OGRCoordinateTransformationH;
OGRSpatialReferenceH OSRNewSpatialReference( const char *);
void   OSRDestroySpatialReference( OGRSpatialReferenceH );
int    OSRReference( OGRSpatialReferenceH );
int    OSRDereference( OGRSpatialReferenceH );
OGRErr OSRImportFromEPSG( OGRSpatialReferenceH, int );
OGRErr OSRImportFromWkt( OGRSpatialReferenceH, char ** );
OGRErr OSRExportToWkt( OGRSpatialReferenceH, char ** );
OGRErr OSRSetAttrValue( OGRSpatialReferenceH hSRS, const char * pszNodePath,
                         const char *pszNewNodeValue );
const char *OSRGetAttrValue( OGRSpatialReferenceH hSRS,
                             const char * pszName, intiChild);
OGRErr OSRSetLinearUnits( OGRSpatialReferenceH, const char *, double );
double OSRGetLinearUnits( OGRSpatialReferenceH, char ** );
int    OSRIsGeographic( OGRSpatialReferenceH );
int    OSRIsProjected( OGRSpatialReferenceH );
int    OSRIsSameGeogCS( OGRSpatialReferenceH, OGRSpatialReferenceH );
int    OSRIsSame( OGRSpatialReferenceH, OGRSpatialReferenceH );
OGRErr OSRSetProjCS( OGRSpatialReferenceH hSRS, const char * pszName );
OGRErr OSRSetWellKnownGeogCS( OGRSpatialReferenceH hSRS,
                               const char * pszName );
OGRErr OSRSetGeogCS( OGRSpatialReferenceH hSRS,
                      const char * pszGeogName,
                      const char *pszDatumName,
                      const char *pszEllipsoidName,
                      double dfSemiMajor,double dfInvFlattening,
                      const char * pszPMName ,
                     double dfPMOffset ,
                      const char * pszUnits,
                      double dfConvertToRadians);
double OSRGetSemiMajor( OGRSpatialReferenceH, OGRErr * );
double OSRGetSemiMinor( OGRSpatialReferenceH, OGRErr * );
double OSRGetInvFlattening( OGRSpatialReferenceH, OGRErr * );
OGRErr OSRSetAuthority( OGRSpatialReferenceH hSRS,
                         const char *pszTargetKey,
                         const char *pszAuthority,
                         int nCode );
OGRErr OSRSetProjParm( OGRSpatialReferenceH, const char *, double );
double OSRGetProjParm( OGRSpatialReferenceH hSRS,
                        const char *pszParmName,
                        double dfDefault,
                        OGRErr * );
OGRErr OSRSetUTM( OGRSpatialReferenceH hSRS, int nZone, int bNorth );
int    OSRGetUTMZone( OGRSpatialReferenceH hSRS, int *pbNorth );
OGRCoordinateTransformationH
OCTNewCoordinateTransformation( OGRSpatialReferenceH hSourceSRS,
                                OGRSpatialReferenceHhTargetSRS );
void OCTDestroyCoordinateTransformation( OGRCoordinateTransformationH );
int OCTTransform(OGRCoordinateTransformationH hCT,
                  int nCount, double *x, double*y, double *z );

2、Python語言接口

class osr.SpatialReference
   def __init__(self,obj=None):
   def ImportFromWkt( self, wkt ):
   def ExportToWkt(self):
   def ImportFromEPSG(self,code):
   def IsGeographic(self):
   def IsProjected(self):
   def GetAttrValue(self, name, child = 0):
   def SetAttrValue(self, name, value):
   def SetWellKnownGeogCS(self, name):
   def SetProjCS(self, name = "unnamed" ):
   def IsSameGeogCS(self, other):
   def IsSame(self, other):
   def SetLinearUnits(self, units_name, to_meters ):
   def SetUTM(self, zone, is_north = 1):

class CoordinateTransformation:
   def __init__(self,source,target):
   def TransformPoint(self, x, y, z = 0):
   def TransformPoints(self, points):

七、內部實現

OGRCoordinateTransformation依賴於PROJ.4庫。所以要使用坐標轉換的內容，GDAL必須在編譯的時候綁定PROJ4才可以用來進行坐標轉換。如果要使用GDAL的坐標轉換，重投影相關的算法，就必須要有PROJ4庫的支持，否則會轉換失敗。關於PROJ4的編譯和GDAL綁定PROJ4的內容，請參考之前的博文。

 

來源： http://blog.csdn.net/jingss_3/article/details/8240328

來自為知筆記(Wiz)