经纬度与UTM坐标的转换

昨天花了一个晚上的时间搞定了经纬度到UTM坐标的转换，UTM到经纬度的转换还没来得及弄，不过有了前者，后者就有底了～
按理说来，这个LL2UTM的转换以前就有中文的网页介绍过，具体网址我记不清了，应该很容易搜到，还有用C语言实现的函数代码，用的是UGSC的方法。我用它来转换经纬度到UTM，画出的地图和实际的很像，就想当然的认为已经搞定了。
昨天因为需要进行逆向转换，也就是UTM2LL，所有在网上搜索，G到了这个网页，作者是Steven Dutch，最后更新是08年4月1号，愚人节哦～里面讲得挺清楚的，有LL2UTM，也有UTM2LL。
LL2UTM有两种方法，一种是美国军方给出的，一种是美国地质测量部给出的，分别简称为Army和USGS。至于UTM2LL，Army的方法牵扯到查表插值，所以只给出了USGS的方法。最让人开心的是Steven还给出了包含转换实现的Excel文件下载！这个东西作用大大的，后面会说到。
我开始没有多想，直接拿我前面转换好的UTM的X、Y坐标填入Excel文件中，转换出来的值却是和原始的经纬度值大不一样！如果UTM2LL的转换没有错的话，那么就是我的LL2UTM出错了！拿原始经纬度值填入Excel文件中转换成UTM坐标，和我转换出来的果然不一样！还差了不少！
仔细看了看Excel文件中转换的中间值，发现我的中央经纬度值就出错了！原先我写的LL2UTM函数需要把中央经纬度值作为输入参数，我手动算118度的中央经纬度居然算出了120，实际上应该是117，改了以后再转换一遍，就差不多了。但看着尚存的误差，我还是有点不甘心。我原先写的函数实际上用的是USGS方法，Steven建议的是用Army的方法，于是我就又用python实现一遍。这次就曲折大了！
按照Steven给出的公式写完函数后，填入原始经纬度值调用，发现出来的UTM值查了好多！仔细检查了好几遍，没错啊！还好Excel文件中有中间过程值显示，就把中间结果都打印出来比较，发现是最后算K*时出的错。看了Excel中的计算公式，发现问题所在了：和网页中的公式不一样！按照Excel中的计算公式改了再试，可以了～～

从经纬度转换为UTM坐标（USGS方法）

def LL2UTM_USGS(a, f, lat, lon, lonOrigin, FN):
'''
** Input：(a, f, lat, lon, lonOrigin, FN)
** a 椭球体长半轴
** f 椭球体扁率 f=(a-b)/a 其中b代表椭球体的短半轴
** lat 经过UTM投影之前的纬度
** lon 经过UTM投影之前的经度
** lonOrigin 中央经度线
** FN 纬度起始点，北半球为0，南半球为10000000.0m
---------------------------------------------
** Output:(UTMNorthing, UTMEasting)
** UTMNorthing 经过UTM投影后的纬度方向的坐标
** UTMEasting 经过UTM投影后的经度方向的坐标
---------------------------------------------
** 功能描述：UTM投影
** 作者： Ace Strong
** 单位： CCA NUAA
** 创建日期：2008年7月19日
** 版本：1.0
** 本程序实现的公式请参考
** "Coordinate Conversions and Transformations including Formulas" p35.
** & http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm
'''

# e表示WGS84第一偏心率,eSquare表示e的平方
eSquare = 2*f - f*f
k0 = 0.9996

# 确保longtitude位于-180.00----179.9之间
lonTemp = (lon+180)-int((lon+180)/360)*360-180
latRad = math.radians(lat)
lonRad = math.radians(lonTemp)
lonOriginRad = math.radians(lonOrigin)
e2Square = (eSquare)/(1-eSquare)

V = a/math.sqrt(1-eSquare*math.sin(latRad)**2)
T = math.tan(latRad)**2
C = e2Square*math.cos(latRad)**2
A = math.cos(latRad)*(lonRad-lonOriginRad)
M = a*((1-eSquare/4-3*eSquare**2/64-5*eSquare**3/256)*latRad
-(3*eSquare/8+3*eSquare**2/32+45*eSquare**3/1024)*math.sin(2*latRad)
+(15*eSquare**2/256+45*eSquare**3/1024)*math.sin(4*latRad)
-(35*eSquare**3/3072)*math.sin(6*latRad))

# x
UTMEasting = k0*V*(A+(1-T+C)*A**3/6
+ (5-18*T+T**2+72*C-58*e2Square)*A**5/120)+ 500000.0
# y
UTMNorthing = k0*(M+V*math.tan(latRad)*(A**2/2+(5-T+9*C+4*C**2)*A**4/24
+(61-58*T+T**2+600*C-330*e2Square)*A**6/720))
# 南半球纬度起点为10000000.0m
UTMNorthing += FN
return (UTMEasting, UTMNorthing)

从经纬度转换到UTM（Army方法）

def LL2UTM_Army(a, b, lat_ll, lon_ll, FN):
'''
** Input：(a, b, lat, lon, FN)
** a 椭球体长半轴
** b 椭球体短半轴
** lat_ll 经过UTM投影之前的纬度(角度为单位)
** lon_ll 经过UTM投影之前的经度(角度为单位)
** FN 纬度起始点，北半球为0，南半球为10000000.0m
---------------------------------------------
** Output:(UTMEasting, UTMNorthing)
** UTMNorthing 经过UTM投影后的纬度方向的坐标
** UTMEasting 经过UTM投影后的经度方向的坐标
---------------------------------------------
** 功能描述：UTM投影
** 作者： Ace Strong
** 单位： CCA NUAA
** 创建日期：2009年1月14日
** 版本：1.0
** 本程序实现的公式请参考
** http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm
'''

lat = math.radians(lat_ll)
lon = math.radians(lon_ll)
lon0_ll = 6*(int(lon_ll/6)+31)-183
k0 = 0.9996
e = math.sqrt(1-b**2/a**2)
e2 = e**2/(1-e**2)
n = (a-b)/(a+b)
rho = a*(1-e**2)/((1-(e*math.sin(lat))**2)**(3.0/2))
nu = a/((1-(e*math.sin(lat))**2)**(1.0/2))
p = (lon_ll-lon0_ll)*3600/10000
sin1 = math.pi/(180*60*60)

A = a*(1 - n + (5.0/4)*(n**2 - n**3) + (81.0/64)*(n**4 - n**5))
B = (3*a*n/2)*(1 - n + (7.0/8)*(n**2 - n**3) + (55.0/64)*(n**4 - n**5))
C = (15*a*n**2/16)*(1 - n + (3.0/4)*(n**2 - n**3))
D = (35*a*n**3/48)*(1 - n + (11.0/16)*(n**2 - n**3))
E = (315*a*n**4/51)*(1 - n)
S = A*lat - B*math.sin(2*lat) + C*math.sin(4*lat) - D*math.sin(6*lat) + E*math.sin(8*lat)

# y
K1 = S*k0
K2 = nu*math.sin(lat)*math.cos(lat)*sin1**2*k0*(100000000)/2
K3 = (sin1**4*nu*math.sin(lat)*math.cos(lat)**3/24)*(5 - math.tan(lat)**2 + 9*e2*math.cos(lat)**2 + 4*e2**2*math.cos(lat)**4)*k0*(10000000000000000)
UTMNorthing = K1 + K2*p**2 + K3*p**4
# 南半球纬度起点为10000000.0m
UTMNorthing += FN
# x
K4 = k0*sin1*nu*math.cos(lat)*10000
K5 = (sin1*math.cos(lat))**3*(nu/6)*(1 - math.tan(lat)**2 + e2*math.cos(lat)**2)*k0*1000000000000
UTMEasting = K4*p + K5*p**3 + 500000
return (UTMEasting, UTMNorthing)