经纬度坐标是一种地理坐标系,记录的经纬度信息,经度范围是(W)-180~180(E),南半球纬度0~90°S 北半球纬度0~90°N。
UTM分带及计算方法,3°分带(高斯克吕格)和6°分带(UTM)
带号计算公式如下:
Zone = int[(Lon - 1.5)/3] + 1 3°分带
Zone = floor(Lon/6)+31 6°分带
中央经线计算公式如下:
L0 = Zone*3 3°分带
L0 = (6*Zone - 3) - 180 6°分带
更多情况可以参考:UTM投影与高斯克吕格投影中分带带号与中央经线经度的计算关系
方法一:直接使用数学方法计算
参考:直接代码搬运进行修改
GPS经纬度坐标和UTM坐标的相互转换
修改代码如下:
from math import cos, sin, tan, sqrt
import math
longitude, latitude=113.4,32.5
def LonLat2UTM(longitude, latitude):
lat = latitude
lon = longitude
kD2R = math.pi / 180.0
ZoneNumber = math.floor(lon / 6.0) + 31
print("ZoneNumber",ZoneNumber)
L0 = (ZoneNumber-31) * 6+3
# L0 = (ZoneNumber-31)*6+6 #中央经线计算
print("中央经线",L0)
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 #东偏移距离 默认500000
N = 0.9996 * N
# E = 0.9996 * E
E = 0.9996 * (E - 500000.0) + 500000
UTME=E
UTMN=N
return UTME,UTMN
print("输入经纬度是:",longitude," ",latitude)
print("输出经纬度是:",LonLat2UTM(longitude, latitude))
输出结果:
输入经纬度是: 113.4 32.5
ZoneNumber 49
中央经线 111
输出经纬度是: (725482.0885903714, 3598397.4113411955)
方法二:使用现有python库计算
使用python库Proj库计算
代码如下:
from pyproj import Proj
# 首先定义要转换的投影坐标系
# p1 = Proj('+proj=utm +zone=38 +datum=WGS84 +units=m +no_defs')
p1 = Proj(proj='utm',zone=49,ellps='WGS84', preserve_units='m')
lon1,lat1 = p1(113.4,32.5) # 将地理坐标转换为投影坐标,地理坐标为WGS84下的坐标
print(lon1,lat1)
输出结果:
725482.0885891931 3598397.411341262
以上两种方法计算得到的结果相同,上面都是采用6°分带计算的坐标,而且false Easting均为-500000.我们将以上结果和ARCGIS中的点进行比较验证。
如果是需要将投影坐标系转换为地理坐标系,可以增加inverse=True来实现。比如:
lon1,lat1 = p1(x0,y0,inverse=True)
其中x0,y0是投影坐标
坐标点转换验证分析
ARCGIS中的坐标结果显示
首先将ARCGIS中的View菜单中Data Frame Properties窗口中的的坐标系统改为UTM进行显示,参数如下:
或者:
显示坐标为:(725,482.0886 m m,3,598,397.4113 m)
总结
通过以上两种方法计算得到的结果和ARCGIS中的结果是一致的。
补充:Albers投影和地理坐标系之间的转换
# Albers等积投影在中国没有已经定义好的,需要自己定义
proj1 = Proj('+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84')
'''
epsg编号通过epsg官网或者arcmap中查询获得,此为最适合中国的albers投影
'''
lon1,lat1 = proj1(98.9406,38.8399) # 将地理坐标转换为投影坐标,地理坐标为WGS84下的坐标
print(lon1,lat1)
lon2,lat2 = proj1(-516480.0192,4184573.034,inverse=True) # 将投影坐标转换为地理坐标
print(lon2,lat2)