[C++知识库] GDAL中的地理坐标系、投影坐标系及其相互转换

开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> C++知识库 -> GDAL中的地理坐标系、投影坐标系及其相互转换 -> 正文阅读

[C++知识库]GDAL中的地理坐标系、投影坐标系及其相互转换

文章目录

- - - - 地理坐标系
        例1：国内常用地理坐标系
        投影坐标系
        例2：国内常用投影坐标系（不推荐使用）
        坐标转换
        例3：地理坐标转投影坐标
        例4：投影坐标转地理坐标

地理坐标系

原理参考这篇文章：
地理坐标系与投影坐标系区别与联系
https://yunxingluoyun.blog.csdn.net/article/details/123970678

例1：国内常用地理坐标系

#include <cstdio>
#include "gdal_priv.h"
#include <iostream>

int main()
{
	OGRSpatialReference oSRS1;
	oSRS1.SetGeogCS("自定义地理坐标系",// 定义地理坐标系名称
		"WGS_1984",// 基准面
		"WGS84 椭球",// 椭球体名称
		SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,// WGS84椭球体的长半轴，WGS84椭球体扁率的倒数
		"Greenwich", 0.0,// 格林尼治子午线
		"degree", 0.0174532925199433);// 角度度量单位
	char* WGS84_WTK = NULL;
	oSRS1.exportToPrettyWkt(&WGS84_WTK);// 以规整的WTK格式输出地理坐标系的信息
	std::cout << WGS84_WTK << std::endl;
	// CGCS2000
	OGRSpatialReference oSRS2;
	oSRS2.SetWellKnownGeogCS("EPSG:4490");
	char* CGCS2000_WTK = NULL;
	oSRS2.exportToPrettyWkt(&CGCS2000_WTK);
	std::cout << CGCS2000_WTK << std::endl;
	// Beijing_1954
	OGRSpatialReference oSRS3;
	oSRS3.SetWellKnownGeogCS("EPSG:4214");
	char* Beijing_1954_WTK = NULL;
	oSRS3.exportToPrettyWkt(&Beijing_1954_WTK);
	std::cout << Beijing_1954_WTK << std::endl;
	// Xian_1980
	OGRSpatialReference oSRS4;
	oSRS4.SetWellKnownGeogCS("EPSG:4610");
	char* Xian_1980_WTK = NULL;
	oSRS4.exportToPrettyWkt(&Xian_1980_WTK);
	std::cout << Xian_1980_WTK << std::endl;

}

结果：

GEOGCS["自定义地理坐标系",
    DATUM["WGS_1984",
        SPHEROID["WGS84 椭球",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST]]
    
GEOGCS["China Geodetic Coordinate System 2000",
    DATUM["China_2000",
        SPHEROID["CGCS2000",6378137,298.257222101,
            AUTHORITY["EPSG","1024"]],
        AUTHORITY["EPSG","1043"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4490"]]
    
GEOGCS["Beijing 1954",
    DATUM["Beijing_1954",
        SPHEROID["Krassowsky 1940",6378245,298.3,
            AUTHORITY["EPSG","7024"]],
        AUTHORITY["EPSG","6214"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4214"]]
   
GEOGCS["Xian 1980",
    DATUM["Xian_1980",
        SPHEROID["IAG 1975",6378140,298.257,
            AUTHORITY["EPSG","7049"]],
        AUTHORITY["EPSG","6610"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4610"]]

投影坐标系

设置横轴墨卡托投影的函数SetTM()有六个参数：

参数	设置
dfCenterLat	一般为0
dfCenterLong	中央经线，决定了是哪一投影带
dfScale	一般为1.0，UTM设置为0.9996
dfFalseEasting	一般为500000，高斯-克吕格前面加上带号
dfFalseNorthing	一般为0

例2：国内常用投影坐标系（不推荐使用）

推荐使用例3中的方式使用坐标系，EPSG请查看网址：

EPSG.io: Coordinate Systems Worldwide
注意：这种方式定义的坐标轴与例3相反，如果使用这样方法，后面的poCT->Transform(1, &x, &y)影改为poCT->Transform(1, &y, &x)。

#include <cstdio>
#include "gdal_priv.h"
#include <iostream>

int main()
{
	OGRSpatialReference oSRS1;
	oSRS1.SetProjCS("UTM 17(WGS84) in northern hemisphere.");// 设置投影坐标系名称
	oSRS1.SetWellKnownGeogCS("EPSG:4326");// 设置地理坐标系
	oSRS1.SetUTM(17, TRUE);//设置UTM投影的投影参数
	char* UTM17_WTK = NULL;
	oSRS1.exportToPrettyWkt(&UTM17_WTK);
	std::cout << UTM17_WTK << std::endl;

	OGRSpatialReference oSRS2;
	oSRS2.SetProjCS("CGCS2000 / Gauss-Kruger CM 117E");
	oSRS2.SetWellKnownGeogCS("EPSG:4490");
	oSRS2.SetTM(0,117,1,500000,0);//起始纬度，中央经线，比例因子，东偏移量，北偏移量
	char* CGCS2000_GK_117E_WTK = NULL;
	oSRS2.exportToPrettyWkt(&CGCS2000_GK_117E_WTK);
	std::cout << "CGCS2000 / Gauss-Kruger CM 117E" << std::endl;
	std::cout <<"是否为投影坐标系：" << oSRS2.IsProjected() << std::endl;
	std::cout << CGCS2000_GK_117E_WTK << std::endl;

	OGRSpatialReference oSRS3;
	oSRS3.SetProjCS("CGCS2000 3 Degree GK Zone 39");
	oSRS3.SetWellKnownGeogCS("EPSG:4490");
	oSRS3.SetTM(0, 117, 1, 39500000, 0);//起始纬度，中央经线，比例因子，东偏移量，北偏移量
	char* CGCS2000_GK_Zone39_WTK = NULL;
	oSRS3.exportToPrettyWkt(&CGCS2000_GK_Zone39_WTK);
	std::cout << "CGCS2000 3 Degree GK Zone 39" << std::endl;
	std::cout << "是否为投影坐标系：" << oSRS3.IsProjected() << std::endl;
	std::cout << CGCS2000_GK_Zone39_WTK << std::endl;

	OGRSpatialReference oSRS4;
	oSRS4.SetProjCS("Beijing 1954 3 Degree GK CM 117E");
	oSRS4.SetWellKnownGeogCS("EPSG:4214");
	oSRS4.SetTM(0, 117, 1, 500000, 0);//起始纬度，中央经线，比例因子，东偏移量，北偏移量
	char* Beijing1954_GK_117E_WTK = NULL;
	oSRS4.exportToPrettyWkt(&Beijing1954_GK_117E_WTK);
	std::cout << "Beijing 1954 3 Degree GK CM 117E" << std::endl;
	std::cout << "是否为投影坐标系：" << oSRS4.IsProjected() << std::endl;
	std::cout << Beijing1954_GK_117E_WTK << std::endl;

	OGRSpatialReference oSRS5;
	oSRS5.SetProjCS("Xian 1980 3 Degree GK Zone 39");
	oSRS5.SetWellKnownGeogCS("EPSG:4610");
	oSRS5.SetTM(0, 117, 1, 39500000, 0);//起始纬度，中央经线，比例因子，东偏移量，北偏移量
	char* Xian_1980_GK_Zone39_WTK = NULL;
	oSRS5.exportToPrettyWkt(&Xian_1980_GK_Zone39_WTK);
	std::cout << "Xian 1980 3 Degree GK Zone 39" << std::endl;
	std::cout << "是否为投影坐标系：" << oSRS5.IsProjected() << std::endl;
	std::cout << Xian_1980_GK_Zone39_WTK << std::endl;
}

结果：

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["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]
    
CGCS2000 / Gauss-Kruger CM 117E
是否为投影坐标系：1
PROJCS["CGCS2000 / Gauss-Kruger CM 117E",
    GEOGCS["China Geodetic Coordinate System 2000",
        DATUM["China_2000",
            SPHEROID["CGCS2000",6378137,298.257222101,
                AUTHORITY["EPSG","1024"]],
            AUTHORITY["EPSG","1043"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4490"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]
    
CGCS2000 3 Degree GK Zone 39
是否为投影坐标系：1
PROJCS["CGCS2000 3 Degree GK Zone 39",
    GEOGCS["China Geodetic Coordinate System 2000",
        DATUM["China_2000",
            SPHEROID["CGCS2000",6378137,298.257222101,
                AUTHORITY["EPSG","1024"]],
            AUTHORITY["EPSG","1043"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4490"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",39500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]
    
Beijing 1954 3 Degree GK CM 117E
是否为投影坐标系：1
PROJCS["Beijing 1954 3 Degree GK CM 117E",
    GEOGCS["Beijing 1954",
        DATUM["Beijing_1954",
            SPHEROID["Krassowsky 1940",6378245,298.3,
                AUTHORITY["EPSG","7024"]],
            AUTHORITY["EPSG","6214"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4214"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]
    
Xian 1980 3 Degree GK Zone 39
是否为投影坐标系：1
PROJCS["Xian 1980 3 Degree GK Zone 39",
    GEOGCS["Xian 1980",
        DATUM["Xian_1980",
            SPHEROID["IAG 1975",6378140,298.257,
                AUTHORITY["EPSG","7049"]],
            AUTHORITY["EPSG","6610"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4610"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",39500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

坐标转换

例3：地理坐标转投影坐标

#include <cstdio>
#include "gdal_priv.h"

int main()
{
    OGRSpatialReference oSourceSRS, oTargetSRS;
    OGRCoordinateTransformation* poCT;
    oSourceSRS.importFromEPSG(4610); // 地理坐标系：Xian 1980
    oTargetSRS.importFromEPSG(2383); // 投影坐标系：Xian 1980 / 3-degree Gauss-Kruger CM 114E
    
    char* Xian_1980 = NULL;
    oSourceSRS.exportToPrettyWkt(&Xian_1980);
    char* Xian_1980_GK_114E_WTK = NULL;
    oTargetSRS.exportToPrettyWkt(&Xian_1980_GK_114E_WTK);
    printf("%s,%s\n", "oSourceSRS:\n", Xian_1980);
    printf("%s,%s\n", "oTargetSRS:\n", Xian_1980_GK_114E_WTK);


    poCT = OGRCreateCoordinateTransformation(&oSourceSRS, &oTargetSRS);
    double x, y;
    x = 38.8; //纬度
    y = 113.6; //经度
    printf("经纬度坐标：%.9lf	%.9lf", x, y);
    if (poCT == NULL || !poCT->Transform(1, &x, &y))
    {
        printf("Transformation failed.\n");
    }
    else
    {
        printf("\n平面坐标：%.9lf	%.9lf", x, y);
    }
    
}

结果：

oSourceSRS:
,GEOGCS["Xian 1980",
    DATUM["Xian_1980",
        SPHEROID["IAG 1975",6378140,298.257,
            AUTHORITY["EPSG","7049"]],
        AUTHORITY["EPSG","6610"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4610"]]
oTargetSRS:
,PROJCS["Xian 1980 / 3-degree Gauss-Kruger CM 114E",
    GEOGCS["Xian 1980",
        DATUM["Xian_1980",
            SPHEROID["IAG 1975",6378140,298.257,
                AUTHORITY["EPSG","7049"]],
            AUTHORITY["EPSG","6610"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4610"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",114],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Northing",NORTH],
    AXIS["Easting",EAST],
    AUTHORITY["EPSG","2383"]]
经纬度坐标：38.800000000        113.600000000
平面坐标：4296379.277209882     465252.023887851

例4：投影坐标转地理坐标

#include <cstdio>
#include "gdal_priv.h"

int main()
{
    OGRSpatialReference oSourceSRS, oTargetSRS;
    OGRCoordinateTransformation* poCT;
    oSourceSRS.importFromEPSG(2383); // Xian 1980 / 3-degree Gauss-Kruger CM 114E
    oTargetSRS.importFromEPSG(4610); // Xian 1980 

    char* Xian_1980_GK_114E_WTK = NULL;
    oSourceSRS.exportToPrettyWkt(&Xian_1980_GK_114E_WTK);
    char* Xian_1980 = NULL;
    oTargetSRS.exportToPrettyWkt(&Xian_1980);

    printf("%s,%s\n", "oSourceSRS:\n", Xian_1980);
    printf("%s,%s\n", "oTargetSRS:\n", Xian_1980_GK_114E_WTK);


    poCT = OGRCreateCoordinateTransformation(&oSourceSRS, &oTargetSRS);
    double x, y;  
    x = 4296379.277209882; //纬度
    y = 465252.023887851; //经度
    printf("\n平面坐标：%.9lf	%.9lf", x, y);
    
    if (poCT == NULL || !poCT->Transform(1, &x, &y))
    {
        printf("Transformation failed.\n");
    }
    else
    {
        printf("\n经纬度坐标：%.9lf  %.9lf", x, y);
    }
    
}

结果：

oSourceSRS:
,GEOGCS["Xian 1980",
    DATUM["Xian_1980",
        SPHEROID["IAG 1975",6378140,298.257,
            AUTHORITY["EPSG","7049"]],
        AUTHORITY["EPSG","6610"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4610"]]
oTargetSRS:
,PROJCS["Xian 1980 / 3-degree Gauss-Kruger CM 114E",
    GEOGCS["Xian 1980",
        DATUM["Xian_1980",
            SPHEROID["IAG 1975",6378140,298.257,
                AUTHORITY["EPSG","7049"]],
            AUTHORITY["EPSG","6610"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4610"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",114],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Northing",NORTH],
    AXIS["Easting",EAST],
    AUTHORITY["EPSG","2383"]]

平面坐标：4296379.277209882     465252.023887851
经纬度坐标：38.800000000  113.600000000