GDAL学习笔记——OGR投影(一)

GDAL学习笔记主要是对官方帮助文档学习的记录,虽然看的有点儿慢,但是还是可以看下去的,希望学习GDAL的同道人静下心来一起学习。学习计划:先学习官方帮助文档中的Tutorials,然后根据实际的应用进一步掌握GDAL的使用。

代码测试环境:
1、Qt 5.6.1(MSVC 2013,64 bit ) + GDAL(MSVC 2013 ,64bit )
2、VS2013(64 bit) + GDAL-2.1.1

定义地理坐标系

在学习这部分内容之前,应该对空间参考(Spatial Reference)、大地水准面基准面(datum)、地图投影(map projection)、地理坐标系(Geographic Coordinate System)和投影坐标系(Projeetion coordinate system)有所了解。

关于坐标系的定义及内容封装在OGRSpatialReference类中,主要的两种坐标系:地理坐标系和投影坐标系。

地理坐标系需要包含基准面、本初子午线角度单位这些信息。下面进行地理坐标系的定义:

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

这里说明下,官网上的最后一个参数是用SRS_UA_DEGREE_CONV宏定义表示的,但是我看了下该函数的原型,最后一个参数是double型的,宏定义如下(表示的是字符串):

#define SRS_UA_DEGREE           "degree"
#define SRS_UA_DEGREE_CONV                  "0.0174532925199433"
#define SRS_UA_RADIAN           "radian"

SetGeogCS函数原型如下:

OGRErr OGRSpatialReference::SetGeogCS(const char *  pszGeogName,
                        const char *    pszDatumName,
                        const char *    pszSpheroidName,
                        double  dfSemiMajor,
                        double  dfInvFlattening,
                        const char *    pszPMName = NULL,
                        double  dfPMOffset = 0.0,
                        const char *    pszAngularUnits = NULL,
                        double  dfConvertToRadians = 0.0 
                        )   

在实际开发中,直接写SRS_UA_DEGREE_CONV是提示错误的(实参与形参不一致),我不知道这部分应该怎么处理?直接写0.0174532925199433吗?希望有知道的大神能告知。

也可以用OGRSpatialReference自带一些标识来表示一些常用的坐标系,比如”NAD27”, “NAD83”, “WGS72” and “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"]]

已知WKT中描述的坐标系,可以利用OGRSpatialReference::importFromWkt()函数来定义OGRSpatialReference对象。

定义投影坐标系

一个投影坐标系(比如UTM,兰伯特等角圆锥投影等)需要以地理坐标系为基础。同时,需要投影方法用于线性位置(米或英尺)与经纬度表示的位置之间的转换。以下代码定义了一个基于WGS84的大地基准椭球体的UTM第17带的投影坐标系统。

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

打印的结果:

PROJCS["UTM 17 (wgs84) in northern hemisphere",
    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"]],
    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],
    UNIT["Meter",1]]

当然OGRSpatialReference不止提供SetUTM()这一种方法,设置横轴墨卡托投影参数可以使用SetTM()函数;设置兰勃特投影参数使用SetLCC()函数;设置墨卡托投影参数使用SetMercator()函数。

上述内容是参考官方的帮助文档,有些地方表达的不准确或不恰当,希望各路英雄好汉批评指正!

©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页