使用GDAL实现坐标转换

在接触这部分内容前，应该对**空间参考（Spatial Reference）、大地水准面、基准面（datum）、地图投影（map projection）、地理坐标系（Geographic Coordinate System）和投影坐标系（Projeetion coordinate system）**有一定了解。OGRSpatialReference类和OGRCoordinateTransformation类主要用来提供定义坐标系统（投影和水准面）和转换坐标。这两个类都基于OpenGIS的坐标转换说明，并且使用Well Known Text格式来进行表述坐标系统。下面来说明如何实现不同坐标系之间的转换：

# **一、定义地理坐标系**

一个地理坐标系统需要包含的信息有一个**大地基准面**（里面含有一个使用长半轴和扁率的倒数来表示的托球体），**一个中央经线**（通常是本初子午线，也就是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”等。

oSRS.SetWellKnownGeogCS("WGS84");

也可以用EPSG数据库中含有的地理坐标系统编码（GSC code）定义坐标系：

oSRS.SetWellKnownGeogCS("EPSG:4326");

为了方便和其他库进行交互，OGRSpatialReference提供了可以和OpenGIS的WKT格式的相互转换的函数。OGRSpatialReference可以使用一个WKT格式文件来进行初始化，也可以将坐标系的信息导出为WKT格式。

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

打印出的结果如下：

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

可以利用OGRSpatialReference::importFromWkt()函数来定义OGRSpatialReference对象。

# **二、定义投影坐标系**

一个投影坐标系统（比如UTM，兰伯特等角圆锥投影等）需要建立在一个地理坐标系统之上，在投影坐标系统中，坐标点使用米或者英尺等长度单位来表示，同时也可以用经纬度的角度坐标来表示。下面将定义一个UTM的第17带的投影坐标系统，基于WGS84的大地基准椭球体。

OGRSpatialReference oPGS;
    oPGS.SetProjCS( "UTM 17 (wgs84) in northern hemisphere" ); // SetProjCS()设置投影坐标系的名称
    oPGS.SetWellKnownGeogCS( "WGS84" ); //SetWellKnownGeogCS()设置地理坐标系
    oPGS.SetUTM( 17, TRUE );    //SetUTM()设置投影变换参数
    oPGS.exportToWkt( &pszWKT);
    printf( "%s\n", pszWKT);

> **说明：**
> （1）首先调用SetProjCS()函数设置投影坐标系统的名称，然后使用函数SetWellKnownGeogCS()指定地理坐标系统，最后调用函数SetUTM()设置投影转换参数信息。完成这些工作之后就定义了一个有效的投影坐标系统。**这里必须要注意定义OGRSpatialReference的顺序！**
> （2）当然OGRSpatialReference不止提供SetUTM()这一种方法，设置横轴墨卡托投影参数可以使用SetTM()函数；设置兰勃特投影参数使用SetLCC()函数；设置墨卡托投影参数使用SetMercator()函数。

# **三、解析坐标系**

当OGRSpatialReference对象被创建时，就可以解析该对象所包含的信息。可以使用OGRSpatialReference提供的IsProjected()和IsGeographic()函数分别判别投影坐标系或地理坐标系是否建立，GetSemiMajor()、GetSemiMinor() 和GetInvFlattening()函数分别获取椭球体的长半轴、短半轴以及扁率的倒数。GetAttrValue()获取PROJCS、GEOGCS、DATUM、SPHEROID和PROJECTION名称的字符串。GetProjParm()函数获取投影的参数信息。GetLinearUnits()函数获取单位类型，并且转换为单位米。

OGRSpatialReference oPGS;
    oPGS.SetWellKnownGeogCS( "WGS84" );
    std::cout<< "SemiMajor->" << oPGS.GetSemiMajor() << std::endl;
    std::cout<< "SemiMinor->" << oPGS.GetSemiMinor() << std::endl;
    std::cout<< "InvFlattening->" << oPGS.GetInvFlattening() << std::endl;

# **三、坐标转换**

OGRCoordinateTransformation类可以用来在不同的坐标系统中进行坐标转换。可以使用函数OGRCreateCoordinateTransformation()创建一个新的坐标转换对象，然后使用OGRCoordinateTransformation::Transform()方法来进行坐标转换。

//Transform()的函数原型
    virtual int OGRCoordinateTransformation::Transform (int nCount,double * x,double * y,double * z = NULL)

对于这四个参数：nCount表示要转换点的个数；x,y,z应该分别是三维坐标点的三个值，z可以为NULL，表示只转换水准面上的点。不考虑高程值。

（1）下面是WGS84与Xian 1980/Gauss-Kruger zone 13之间转换的测试代码：

OGRSpatialReference oSourceSRS, oTargetSRS;
    OGRCoordinateTransformation *poCT;
    double x, y;
    oSourceSRS.importFromEPSG(4326);
    WGS84oTargetSRS.importFromEPSG(2327);//Xian 1980/Gauss-Kruger zone 13
    poCT = OGRCreateCoordinateTransformation(&oSourceSRS,&oTargetSRS );
    x = 10.0;
    y = 10.0;
    if (poCT == NULL || !poCT->Transform( 1, &x, &y ))
        cout << "Transformation failed.\n" << endl;
    elsep
        cout << "(%f,%f) -> (%f,%f)\n",10.0,10.0,x, y <<endl;

运行结果： 
(10.000000,10.000000) -> (4381699.110753,2447192.109090)

关于代码中importFromEPSG函数的EPSG代号，可以点击这里查看，或者点我呀也可以查看。

（2）该函数还可以对三维点进行转换，根据不同的椭圆球及基准面自动调整高程值。如果没有Z值，OGR则假设转换的点都是在水准面上。接下来的代码演示了地理坐标系与投影坐标系之间的转换。

OGRSpatialReference    oUTM, *poLatLong;OGRCoordinateTransformation *poTransform;oUTM.SetProjCS( "UTM 17 / WGS84" );oUTM.SetWellKnownGeogCS( "WGS84" );oUTM.SetUTM( 17 );poLatLong = oUTM.CloneGeogCS();poTransform = OGRCreateCoordinateTransformation( &oUTM, poLatLong );double x, y, z;x = 0.0;y = 0.0;z = 50.0;if ( poTransform == NULL || !poTransform->Transform( 1, &x, &y, &z ) )printf( "Transformation failed.\n" );else{printf( "(%f,%f,%f) -> (%f,%f,%f)\n",0.0,0.0,50.0,x, y, z );}
运行结果： 
(0.000000,0.000000,50.000000) -> (-85.488744,0.000000,50.000000)

## **1. 什么是C++**1、

An's Blog

Navigation

Recent Posts

Friend Links