Proj4 坐标转换（以北京54坐标转WGS84投影坐标系为例）

最新推荐文章于 2024-05-28 21:13:10 发布

去往火星

最新推荐文章于 2024-05-28 21:13:10 发布

阅读量215

点赞数 1

分类专栏：遥感开发学习文章标签： c++

本文链接：https://blog.csdn.net/qq_35662333/article/details/137635702

版权

遥感开发学习专栏收录该内容

60 篇文章 2 订阅

订阅专栏

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