您的位置:首页 > 编程语言 > C语言/C++

用Proj.4进行坐标系转换(以北京54坐标系转WGS84投影坐标系为例)

2016-06-17 18:02 645 查看
示例 代码如下:
// convertor.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include "convertor.h"
#include <projects.h>

#ifdef _DEBUG
#define new DEBUG_NEW
#endif

// 唯一的应用程序对象

void output_coordinates(const double coords[3][3], const char* msg, bool degtorad = false);

CWinApp theApp;

using namespace std;

int _tmain(int argc, TCHAR* argv[], TCHAR* envp[])
{
int nRetCode = 0;

HMODULE hModule = ::GetModuleHandle(NULL);

if (hModule != NULL)
{
if (!AfxWinInit(hModule, NULL, ::GetCommandLine(), 0))
{
_tprintf(_T("错误:  MFC 初始化失败\n"));
nRetCode = 1;
}
else
{
}
}
else
{
_tprintf(_T("错误:  GetModuleHandle 失败\n"));
nRetCode = 1;
}

projPJ lcc = pj_init_plus(" +proj=lcc +x_0=0 +y_0=0 +lat_0=0 +lon_0=105 +lat_1=30 +lat_2=62 +a=6378245 +b=6356863.01877305 +no_defs");
//projPJ lcc = pj_init_plus(" +proj=tmerc +lat_0=0 +lon_0=105 +k=1 +x_0=500000 +y_0=0 +a=6378245 +ellps=krass +towgs84=15.8,-154.4,-82.3,0,0,0,0 +units=m +no_defs");
projPJ lonlat1 = pj_init_plus(" +proj=longlat +ellps=krass +towgs84=15.8,-154.4,-82.3,0,0,0,0 +no_defs");
projPJ lonlat2 = pj_init_plus(" +proj=longlat +datum=WGS84 +no_defs");
projPJ merc = pj_init_plus(" +proj=merc +x_0=0 +y_0=0 +lon_0=0 +lat_1=0 +datum=WGS84 +no_defs");
double xyz1[3][3] = { { 1000000, 1100000, 1 }, { 2000000, 2200000, 2 }, { 3000000, 3300000, 3 } };
double xyz2[3][3] = { { 1000000, 1100000, 1 }, { 2000000, 2200000, 2 }, { 3000000, 3300000, 3 } };
int ret = 0;

if (lcc == nullptr || lonlat1 == nullptr || lonlat2 == nullptr || merc == nullptr)
{
printf_s("坐标系初始化失败。\n");
goto zanting;
}

output_coordinates(xyz1, "北京54投影坐标系");

for (int i = 0; i < 3; i++)
ret = pj_transform(lcc, lonlat1, 1, 1, &xyz1[i][0], &xyz1[i][1], &xyz1[i][2]);
output_coordinates(xyz1, "北京54投影坐标系转大地坐标系", true);
output_coordinates(xyz1, "北京54投影坐标系转大地坐标系");

for (int i = 0; i < 3; i++)
pj_geodetic_to_geocentric(lcc->a, lcc->es, 1, 1, &xyz1[i][0], &xyz1[i][1], &xyz1[i][2]);
output_coordinates(xyz1, "大地坐标系转空间直角坐标系");

for (int i = 0; i < 3; i++)
pj_geocentric_to_geodetic(merc->a, merc->es, 1, 1, &xyz1[i][0], &xyz1[i][1], &xyz1[i][2]);
output_coordinates(xyz1, "空间直角坐标系转大地坐标系", true);

for (int i = 0; i < 3; i++)
{
xyz1[i][0] *= DEG_TO_RAD;
xyz1[i][1] *= DEG_TO_RAD;
ret = pj_transform(lonlat2, merc, 1, 1, &xyz1[i][0], &xyz1[i][1], &xyz1[i][2]);
}
output_coordinates(xyz1, "大地坐标系转Mecator投影坐标系");
for (int i = 0; i < 3; i++)
ret = pj_transform(lcc, merc, 1, 1, &xyz2[i][0], &xyz2[i][1], &xyz2[i][2]);
output_coordinates(xyz2, "北京54投影坐标系直接转Mecator投影坐标系");
zanting:
_tsystem(_T("pause"));
return nRetCode;
}

void output_coordinates(const double coords[3][3], const char* msg, bool degtorad)
{
printf_s("=================================================\n");
printf_s(msg);
printf_s("\n");
for (int i = 0; i < 3; i++)
{
if (degtorad)
printf_s("%16f,%16f,%16f\n", coords[i][0] / DEG_TO_RAD, coords[i][1] / DEG_TO_RAD, coords[i][2]);
else
printf_s("%16f,%16f,%16f\n", coords[i][0], coords[i][1], coords[i][2]);
}
printf_s("=================================================\n");
}
代码中有很多中文标注,应该不用多讲吧。
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签:  C++ GIS Proj GDAL