用Proj.4进坐标系转换(不涉及椭球转换)


ProjNet精度有限,支持的坐标系也有限,用WKT表示的坐标系转换参数远没有proj4表示的清爽和一目了然。于是试试proj.4。

proj.4本身不提供坐标系转换参数,椭球转换参数只保存了几种,在pj_data.c代码文件中。需要用7参数指定椭球转换参数的,可以增加+towgs84参数,7个参数的数值用半角逗号隔开就行了。

在转换前调用

projPJ pj_init_plus(const char* proj4text)


以实始化坐标系,分别初始化源坐标系和目标坐标系,

然后调用:

int pj_transform(projPJ from,projPJ to,int count,int offset,double *x,double *y,double *z)

进行坐标转换。下面是一个示例:

// 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");
}

本示例使用源坐标系是D_Beijing_1954坐标系,中央经线是105度。目标坐标系是WGS84。代码有很多中文标注。应该不用多讲吧。

智能推荐

注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
© 2014-2019 ITdaan.com 粤ICP备14056181号  

赞助商广告