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");
}
本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。