GPS坐标转换：经纬度转UTM坐标（Matlab程序）

GPS的输出一般为经纬度坐标，如：北纬（N）41度30分5秒， 东经（E）114度10分30秒。
UTM（Universal Transverse Mercator）坐标是一个投影坐标系，将地球分为60个区，单位为m。
二者的转换可以利用商业软件，如Global Mapper等。下面给出了经纬度转UTM的matlab程序。

    转换公式：（只针对北半球，中国地区计算方法）
   

    Matlab 程序：
   

    运行结果：
    单位为km
   

Global Mapper结果对比：

结果与商业软件在毫米位上保持一致，验证程序的有效性。
---------------------  
作者：xyJupiter  
来源：CSDN  
原文：https://blog.csdn.net/Action_MYX/article/details/83421594  
版权声明：本文为博主原创文章，转载请附上博文链接！

lat = 40.0691643333333;
lon = 116.242161333333;  

kD2R = pi/180;
ZoneNumber = floor((lon - 1.5) / 3.0) + 1;
  L0 = ZoneNumber * 3.0;

  a = 6378137.0;
  F = 298.257223563;
  f = 1 / F;
  b = a * (1 - f);
  ee = (a * a - b * b) / (a * a);
  e2 = (a * a - b * b) / (b * b);
  n = (a - b) / (a + b); n2 = (n * n); n3 = (n2 * n); n4 = (n2 * n2); n5 = (n4 * n);
  al = (a + b) * (1 + n2 / 4 + n4 / 64) / 2.0;
  bt = -3 * n / 2 + 9 * n3 / 16 - 3 * n5 / 32.0;
  gm = 15 * n2 / 16 - 15 * n4 / 32;
  dt = -35 * n3 / 48 + 105 * n5 / 256;
  ep = 315 * n4 / 512;
  B = lat * kD2R;
  L = lon * kD2R;
  L0 = L0 * kD2R;
  l = L - L0;cl = (cos(B) * l);cl2 = (cl * cl); cl3 = (cl2 * cl); cl4 = (cl2 * cl2); cl5 = (cl4 * cl); cl6 = (cl5 * cl); cl7 = (cl6 * cl); cl8 = (cl4 * cl4);
  lB = al * (B + bt * sin(2 * B) + gm * sin(4 * B) + dt * sin(6 * B) + ep * sin(8 * B));
  t = tan(B); t2 = (t * t); t4 = (t2 * t2); t6 = (t4 * t2);
  Nn = a / sqrt(1 - ee * sin(B) * sin(B));
  yt = e2 * cos(B) * cos(B);
  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;
  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

E=4.353761057229367e+05;

N=4.435708949468550e+06;