GPS经纬度坐标转UTM坐标(c++)

UTM的简单介绍
我们在地理课上学过用经纬度表示在地球上的坐标，这种表示方法适合在一个三维的地球曲面上
描述一个坐标，但是如果我们想在二维的平面上描述一个地理坐标应该怎么做呢？

1. 首先是将地球曲面投影到平面上
如下图在地球中心为投影中心，向圆柱形进行射线的投影，然后将圆柱形的侧面进行展开，就可以得到一个平面投影的坐标。

如下图所示就是一幅世界的平面地图，正如我们平时看到的地图

2.将平面分成许多个cell，叫做zone
如下图，分成了n个网格，以向东的方向为x轴，向北的方向为y轴，建立坐标系，每个小格子的中间的x坐标值都是500,000m

因为地球的是圆的，因此每个格子x轴方向的长度是不一样的，如下图

北方向是从南极指向北极，以赤道为分界线，坐标变化为0mS-1,000,000mS，再到0mN-1,000,000mN。其中S和N分别代表南半球和北半球，单位为:m

UTM的东北坐标只有两个数字，还需要zone的区间号才能转换为经纬度坐标因此准确的UTM坐标包括三个参数，区间+东北坐标。如下图

从经纬度转换为UTM坐标的E,N，然后zone编号就得到了唯一的UTM坐标了

经纬度转UTM的东北坐标如下

void LonLat2UTM(double longitude, double latitude, double& UTME, double& UTMN)
{
	double lat = latitude;
	double lon = longitude;

	double kD2R = PI / 180.0;
	double ZoneNumber = floor((lon - 1.5) / 3.0) + 1;
	double L0 = ZoneNumber * 3.0;

	double a = 6378137.0;
	double F = 298.257223563;
	double f = 1 / F;
	double b = a * (1 - f);
	double ee = (a * a - b * b) / (a * a);
	double e2 = (a * a - b * b) / (b * b);
	double n = (a - b) / (a + b); 
	double n2 = (n * n); 
	double n3 = (n2 * n); 
	double n4 = (n2 * n2); 
	double n5 = (n4 * n);
	double al = (a + b) * (1 + n2 / 4 + n4 / 64) / 2.0;
	double bt = -3 * n / 2 + 9 * n3 / 16 - 3 * n5 / 32.0;
	double gm = 15 * n2 / 16 - 15 * n4 / 32;
	double dt = -35 * n3 / 48 + 105 * n5 / 256;
	double ep = 315 * n4 / 512;
	double B = lat * kD2R;
	double L = lon * kD2R;
	L0 = L0 * kD2R;
	double l = L - L0; 
	double cl = (cos(B) * l); 
	double cl2 = (cl * cl); 
	double cl3 = (cl2 * cl); 
	double cl4 = (cl2 * cl2); 
	double cl5 = (cl4 * cl); 
	double cl6 = (cl5 * cl); 
	double cl7 = (cl6 * cl); 
	double cl8 = (cl4 * cl4);
	double lB = al * (B + bt * sin(2 * B) + gm * sin(4 * B) + dt * sin(6 * B) + ep * sin(8 * B));
	double t = tan(B); 
	double t2 = (t * t); 
	double t4 = (t2 * t2); 
	double t6 = (t4 * t2);
	double Nn = a / sqrt(1 - ee * sin(B) * sin(B));
	double yt = e2 * cos(B) * cos(B);
	double N = lB;
	N = N + t * Nn * cl2 / 2;
	N = N + t * Nn * cl4 * (5 - t2 + 9 * yt + 4 * yt * yt) / 24;
	N = N + t * Nn * cl6 * (61 - 58 * t2 + t4 + 270 * yt - 330 * t2 * yt) / 720;
	N = N + t * Nn * cl8 * (1385 - 3111 * t2 + 543 * t4 - t6) / 40320;
	double E = Nn * cl;
	E = E + Nn * cl3 * (1 - t2 + yt) / 6;
	E = E + Nn * cl5 * (5 - 18 * t2 + t4 + 14 * yt - 58 * t2 * yt) / 120;
	E = E + Nn * cl7 * (61 - 479 * t2 + 179 * t4 - t6) / 5040;
	E = E + 500000;
	N = 0.9996 * N;
	E = 0.9996 * (E - 500000.0) + 500000.0;

	UTME = E;
	UTMN = N;
}

转换结果可以利用软件：UTM Coordinate Converter进行验证

1.https://www.youtube.com/watch?v=LcVlx4Gur7I
2.https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system