使用Proj库进行大地坐标转空间坐标投影坐标的一个示例

Posted 2020-08-22 尘世俗人

最近研究了proj库的使用，自己写了一个小demo，仅供参考。

void demoPROJ()
{
	const char* wgs84 = "+proj=tmerc +ellps=WGS84 +lon_0=117e +x_0=500000 +y_0=0 +k=1.0";//+datum=WGS84 

	projPJ pj;

	if (!(pj = pj_init_plus(wgs84)))
	{
		exit(-1);
	}

	projUV pt[2] = {{116.987, 39.2333}, {116.987654321, 39.453}};

	for (int i = 0; i < 2; i++)
	{
		cout << "第" << i + 1 << "个点的转换：\n";
		projUV	pjPt;
		pt[i].u *= DEG_TO_RAD;
		pt[i].v *= DEG_TO_RAD;
		pjPt = pj_fwd(pt[i], pj);
		cout.setf(ios::fixed);
		cout.precision(4);
		cout.width(12);
		cout << "\t经纬度转换为坐标\n\t";
		cout << pjPt.u << "\t" << pjPt.v << endl;

		pjPt = pj_inv(pjPt, pj);
		cout.setf(ios::fixed);
		cout.precision(9);
		cout.width(12);
		cout << "\t坐标转换为经纬度\n\t";
		cout << pjPt.u * RAD_TO_DEG << "\t" << pjPt.v * RAD_TO_DEG << endl;
	}
	cout << endl;

	double wgs84_a = 6378137.0;
	double wgs84_e = 0.00669437999013; // 在这里e是第一偏心率的平方

	projUVW pt3[2] = {{116.5164884833, 39.7663036028,  23.220}, {116.987654321, 39.453, 18}};
	for (int i = 0; i < 2; i++)
	{
		cout << "第" << i + 3 << "个点的转换：\n";
		pt3[i].u *= DEG_TO_RAD;
		pt3[i].v *= DEG_TO_RAD;

		pj_geodetic_to_geocentric(wgs84_a, wgs84_e, 1, 1, &pt3[i].u, &pt3[i].v, &pt3[i].w);
		cout.setf(ios::fixed);
		cout.precision(4);
		cout.width(12);
		cout << "\t经纬度转换为空间坐标\n\t";
		cout << pt3[i].u << "\t" << pt3[i].v << "\t" << pt3[i].w << endl;

		pj_geocentric_to_geodetic(wgs84_a, wgs84_e, 1, 1, &pt3[i].u, &pt3[i].v, &pt3[i].w);
		cout.setf(ios::fixed);
		cout.precision(9);
		cout.width(12);
		cout << "\t空间坐标转换为经纬度\n\t";
		cout << pt3[i].u * RAD_TO_DEG << "\t" << pt3[i].v * RAD_TO_DEG << "\t" ;
		cout.precision(4);
		cout << pt3[i].w << endl;
	}

	pj_free(pj);
}