地图坐标系是用于描述地图上位置的数学模型。它可以用来表示地球表面上的任意一个点，使得这个点的位置可以在地图上精确定位。不同的地图坐标系采用不同的基准面和投影方式，因此会有不同的坐标系参数，不同的坐标系之间也需要进行坐标转换。在地图制图、导航、GIS等领域中，地图坐标系是一个非常重要的概念。
本文介绍了地图坐标系的定义、分类、投影方式、基准面、坐标表示和转换方法等基本概念和原理，以及常见的地图坐标系的特点和应用场景。本文还给出了一些示例代码，以展示如何使用不同的地图API进行坐标系转换。
## 地球经纬度坐标系
地球经纬度坐标系是一种球面坐标系，通常用于描述地球表面上的位置。它将地球分为两个半球，分别为东半球和西半球，以及南半球和北半球，然后通过经度和纬度来确定在球面上的位置。
经度是指在地球上从南极到北极的虚拟线，也称为子午线。经度的取值范围是-180到180度，其中0度经线是本初子午线，位于英国伦敦的皇家格林威治天文台内。东经表示为正数，西经表示为负数。
纬度是指从地球中心到地球表面的一条线，垂直于经度线。纬度的取值范围是-90到90度，其中赤道是0度纬线，南纬表示为负数，北纬表示为正数。
地球经纬度坐标系的优点是易于理解和计算，但由于地球是一个球体，所以在地图制图和空间分析中使用时需要进行坐标变换。
### 示例代码
以下是使用百度地图JS API进行经纬度与容器像素坐标互换的示例代码²：
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<title>容器像素与经纬度互换</title>
<!-- 引入百度地图JS API -->
<script type="text/javascript" src="http://api.map.baidu.com/api?v=3.0&ak=您的密钥"></script>
</head>
<body>
<!-- 创建一个容器用于显示地图 -->
<div id="container" style="width: 600px; height: 400px;"></div>
<!-- 创建一个按钮用于触发事件 -->
<button id="btn">获取当前鼠标位置</button>
<!-- 创建一个文本框用于显示结果 -->
<input id="result" type="text" style="width: 500px;">
<script type="text/javascript">
// 初始化地图对象
var map = new BMap.Map("container");
// 设置中心点和缩放级别
var point = new BMap.Point(116.404, 39.915);
map.centerAndZoom(point, 15);
// 启用滚轮放大缩小
map.enableScrollWheelZoom(true);
// 获取按钮对象
var btn = document.getElementById("btn");
// 获取文本框对象
var result = document.getElementById("result");
// 给按钮绑定点击事件
btn.onclick = function () {
// 获取当前鼠标在地图容器中的像素坐标
var pixel = map.getDefaultCursor();
// 将像素坐标转换为经纬度坐标
var lnglat = map.pixelToPoint(pixel);
// 将经纬度坐标转换为像素坐标
var pixel2 = map.pointToPixel(lnglat);
// 将结果显示在文本框中
result.value = "鼠标在地图容器中的像素坐标为：" + pixel.x + "," + pixel.y + "\n" +
"转换后的经纬度坐标为：" + lnglat.lng + "," + lnglat.lat + "\n" +
"再次转换后的像素坐标为：" + pixel2.x + "," + pixel2.y;
};
</script>
</body>
</html>## 平面直角坐标系（如UTM坐标系）
平面直角坐标系（Planar Coordinate Systems）是一种基于平面直角坐标系的坐标系统，通常用于相对较小的地理区域，例如城市、州、国家等范围内。在平面直角坐标系中，地球表面被划分为网格，每个位置都有唯一的坐标值。
UTM坐标系（Universal Transverse Mercator）是平面直角坐标系中最常用的坐标系统之一。UTM坐标系把地球表面划分成60个纵向带和6个横向带，每个带对应一个投影系统。其中纵向带覆盖从经度0度到经度360度，横向带则依据赤道线和两极线的位置进行划分。每个带内的坐标系统采用横向墨卡托投影或是纵向墨卡托投影，将地球表面投影到平面上，并使用假东、假北坐标表示位置。
相对于经纬度坐标系，平面直角坐标系在短距离内能够提供更高的精度，通常用于测量、地图绘制等领域。但是，由于平面直角坐标系在较长距离内会出现较大的投影误差，因此并不适用于大范围的测量和导航等应用。
### 示例代码
以下是使用高德地图JS API进行经纬度与UTM平面直角坐标互换的示例代码：
<!doctype html>
<html>
<head>
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="initial-scale=1.0, user-scalable=no, width=device-width">
<title>经纬度与平面地图像素互换</title>
<!-- 引入高德地图JS API -->
<script type="text/javascript" src="https://webapi.amap.com/maps?v=1.4.15&key=您的key值"></script>
</head>
<body>
<!-- 创建一个容器用于显示地图 -->
<div id="container" style="width: 600px; height: 400px;"></div>
<!-- 创建一个按钮用于触发事件 -->
<button id="btn">获取当前鼠标位置</button>
<!-- 创建一个文本框用于显示结果 -->
<input id="result" type="text" style="width: 500px;">
<script type="text/javascript">
// 初始化地图对象
var map = new AMap.Map("container", {
zoom: 3 // 设置缩放级别
});
// 获取按钮对象
var btn = document.getElementById("btn");
// 获取文本框对象
var result = document.getElementById("result");
// 给按钮绑定点击事件
btn.onclick = function () {
// 获取当前鼠标在地图容器中的像素坐标
var pixel = map.getDefaultCursor();
// 将像素坐标转换为经纬度坐标
var lnglat = map.containerToLngLat(pixel);
// 将经纬度坐标转换为UTM平面直角坐标
var utm = AMap.LngLat.UTM(lnglat);
// 将UTM平面直角坐标转换为经纬度坐标
var lnglat2 = AMap.LngLat.UTMToLngLat(utm);
// 将结果显示在文本框中
result.value = "鼠标在地图容器中的像素坐标为：" + pixel.x + "," + pixel.y + "\n" +
"转换后的经纬度坐标为：" + lnglat.lng + "," + lnglat.lat + "\n" +
"转换后的UTM平面直角坐标为：" + utm.x + "," + utm.y + "\n" +
"再次转换后的经纬度坐标为：" + lnglat2.lng + "," + lnglat2.lat;
};
</script>
</body>
</html>## 地方坐标系（如高斯-克吕格坐标系）
地方坐标系（Local Coordinate System）是在地球表面上某个点处定义的坐标系，其坐标轴一般是平面直角坐标系，由局部大地基准确定。它与全球坐标系的转换需要考虑地球的椭球形状和大地基准的差异。
其中，高斯-克吕格坐标系（Gauss-Krüger Coordinate System），也叫高斯投影坐标系，是常见的地方坐标系之一。在高斯-克吕格坐标系中，地球被分成若干个带状区域，每个区域内可以通过一个平面直角坐标系来表示坐标。在中国，高斯-克吕格投影一般采用六度带，每个带宽度为 6 度，以中央经线为中心，向东西两边各三度。
由于高斯-克吕格坐标系是地方坐标系，因此其坐标原点和带宽度的选取是与实际地理情况相关的。在中国，高斯-克吕格坐标系一般采用 1954 年版北京坐标系（BJS54）作为大地基准，因此需要将 WGS-84 坐标系下的经纬度坐标先转换到 BJS54 坐标系下，然后再将其转换到高斯-克吕格坐标系下。
与经纬度坐标系相比，高斯-克吕格坐标系具有以下优点：
- 坐标值易于处理和计算，且可表示为实数形式；
- 在较小的区域内，投影误差较小，精度较高；
- 各地图投影的形状、大小、比例等参数是已知的，具有较好的相对精度和可比性。
不过，由于高斯-克吕格坐标系是局部的，其不适用于跨越多个投影带的情况。此外，在大范围地图上使用高斯-克吕格投影会引入较大的投影误差，因此在较大区域内使用 WGS-84 坐标系更为合适。
### 示例代码
以下是使用Python进行WGS-84与高斯-克吕格坐标系互换的示例代码：
# 导入math模块
import math
# 定义WGS-84椭球体参数
a = 6378137.0 # 长半轴
b = 6356752.3142 # 短半轴
f = (a - b) / a # 扁率
e2 = f * (2 - f) # 第一偏心率平方
e12 = e2 / (1 - e2) # 第二偏心率平方
# 定义BJS54大地基准参数
dx = -31.4 # X轴平移量（米）
dy = 144.3 # Y轴平移量（米）
dz = -74.8 # Z轴平移量（米）
rx = 1.735 / 3600 * math.pi / 180 # X轴旋转角（弧度）
ry = -3.041 / 3600 * math.pi / 180 # Y轴旋转角（弧度）
rz = -1.489 / 3600 * math.pi / 180 # Z轴旋转角（弧度）
m = 0.00178 # 尺度因子
# 定义高斯投影中所选用的投影带，如北京为39带
ProjNo = 39
# 定义高斯投影的一些常量
ZoneWide = 6 # 带宽度（度）
PI = math.pi # 圆周率
Rad = PI / 180 # 弧度单位
# 定义一个函数，将WGS-84经纬度坐标转换为BJS54经纬度坐标
def WGS84ToBJS54(lon, lat):
# 将经纬度转换为弧度
lon_rad = lon * Rad
lat_rad = lat * Rad
# 将经纬度转换为XYZ坐标
N = a / math.sqrt(1 - e2 * math.sin(lat_rad) ** 2)
X = N * math.cos(lat_rad) * math.cos(lon_rad)
Y = N * math.cos(lat_rad) * math.sin(lon_rad)
Z = N * (1 - e2) * math.sin(lat_rad)
# 将XYZ坐标进行七参数转换，得到BJS54下的XYZ坐标
X2 = dx + (1 + m) * (X + rz * Y - ry * Z)
Y2 = dy + (1 + m) * (-rz * X + Y + rx * Z)
Z2 = dz + (1 + m) * (ry * X - rx * Y + Z)
# 将BJS54下的XYZ坐标转换为经纬度坐标
lon2_rad = math.atan(Y2 / X2)
Bf = math.atan(Z2 / ((1 - e2) * math.sqrt(X2 ** 2 + Y2 ** 2)))
Nf = a / math.sqrt(1 - e2 * math.sin(Bf) ** 2)
lat2_rad = math.atan((Z2 + Nf * e2 * math.sin(Bf)) / math.sqrt(X2 ** 2 + Y2 ** 2))
lon2_deg = lon2_rad / Rad
lat2_deg = lat2_rad / Rad
return lon2_deg, lat2_deg
# 定义一个函数，将BJS54经纬度坐标转换为高斯-克吕格平面直角坐标
def BJS54ToGauss(lon, lat):
# 将经纬度转换为弧度
lon_rad = lon * Rad
lat_rad = lat * Rad
# 计算中央子午线经度和带号
ZoneWideRad = ZoneWide * Rad
ProjNoRad = ProjNo * ZoneWideRad
ProjNoDeg = ProjNo * ZoneWide
Lon0Rad = ProjNoRad - ZoneWideRad / 2
Lon0Deg = ProjNoDeg - ZoneWide / 2
Lon0SecDeg = Lon0Deg * 3600
# 计算经纬差
e = lon_rad - Lon0Rad
# 计算子午线弧长
a0 = 1 + (3 / 4) * e2 + (45 / 64) * e2 ** 2 + (175 / 256) * e2 ** 3 + (11025 / 16384) * e2 ** 4 + (43659 / 65536) * e2 ** 5 + (693693 / 1048576) * e2 ** 6 + (19324305 / 29360128) * e2 ** 7 + (4927697775 / 7516192768) * e2 ** 8
a2 = -(3 / 4) * e2 - (15 / 16) * e2 ** 2 - (525 / 512) * e2 ** 3 - (2205 / 2048) * e2 ** 4 - (72765 / 65536) * e2 ** 5 - (297297 / 262144) * e2 ** 6 - (135270135 / 117440512) * e2 ** 7 - (547521975 / 469762048) * e2 ** 8
a4 = (15 / 64) * e2 ** 2 + (105 / 256) * e2 ** 3 + (2205 / 4096) * e2 ** 4 + (10395 / 16384) * e2 ** 5 + (1486485 / 2097152) * e2 ** 6 + (45090045 / 58720256) * e2 ** 7 + (766530765 / 939524096) * e2 ** 8
a6 = -(35 / 512) * e2 ** 3 - (315 / 2048) * e2 ** 4 - (31185 / 131072) * e2 ** 5 - (165165 / 524288) * e2 ** 6 - (45090045 / 117440512) * e2 **7 - (209053845 /335544320) *e2**8
a8 = (315/16384)*e2**4+(3465/65536)*e2**5+(99099/1048576)*e2**6+(4099095/29360128)*e2**7+(348423075/1879048192)*e2**8
a10 = -(693/131072)*e2**5-(17325/262144)*e2**6-(99099/524288)*e2**7-(4099095/117440512)*e2**8-(26801775/469762048)*e2**9
a12 = (3003/2097152)*e2**6+(315315/58720256)*e2**7+(11486475/939524096)*e2**8+(90810725/1879048192)*e2**9
X = a0*a*lat_rad+a1*a*math.sin(1*lat_rad)*math.cos(1*lat_rad)+a3*a*math.sin(3*lat_rad)*math.cos(3*lat_rad)+a5*a*math.sin(5*lat_rad)*math.cos(5*lat_rad)+a7*a*math.sin(7*lat_rad)*math.cos(7*lat_rad)+a9*a*math.sin(9*lat_rad)*math.cos(9*lat_rad)+a11*a*math.sin(11*lat_rad)*math.cos(11*lat_rad)
# 计算高斯平面坐标
N = a/math.sqrt(1-e*e*math.sin(lat_rad)**)
t = math.tan(lat_rad)
C = ee*math.cos(lat_rad)**)
A = math.cos(lat_rad)*e
m0 = Lon0SecDeg
m1 = Lon0SecDeg+300000
m = m0+ProjNo*1000000+500000
y = X+N*t*(A*A/1!+A*A*A*A/3!+A*A*A*A*A*A/5!+A*A*A*A*A*A*A*A/7!+A*A*A*A*A*A*A*A*A*A/9!+A*A*A*A*A*A*A*A*A*A*A*A/11!+A*A*A*A*A*A*A*A*A*A*A*A*A*/13!)
x = N*(A+A*A*A/6!+A*A*A*A*A/24!+A*A*A*A*A*A*A/120!+A*A*A*A*A*A*A*A*A/720!+A*A*A*A*A*A*A*A*A*A*/5040!+A*A*A*A*A*A*A*A*A*A/40320!+A*A*A*A*A*A*A*A*A*A*A/362880!+A*A*A*A*A*A*A*A*A*A*A*A*/3628800!)*t+C*(A*A/2!+A*A*A*A/12!+A*A*A*A*A*A/360!+A*A*A*A*A*A*A*/20160!+A*A*A*A*A*A* A* A/151200!+A* A* A* A* A* A* A* A* A/1814400!+A* A* A* A* A* A* A* A* A* A/23950080!)*math.sin(lat_rad)*math.cos(lat_rad)
# 返回高斯平面坐标
return x, y
# 定义一个函数，将高斯-克吕格平面直角坐标转换为BJS54经纬度坐标
def GaussToBJS54(x, y):
# 计算中央子午线经度和带号
ZoneWideRad = ZoneWide * Rad
ProjNoRad = ProjNo * ZoneWideRad
ProjNoDeg = ProjNo * ZoneWide
Lon0Rad = ProjNoRad - ZoneWideRad / 2
Lon0Deg = ProjNoDeg - ZoneWide / 2
Lon0SecDeg = Lon0Deg * 3600
# 计算子午线弧长
a0 = 1 + (3 / 4) * e2 + (45 / 64) * e2 ** 2 + (175 / 256) * e2 ** 3 + (11025 / 16384) * e2 ** 4 + (43659 / 65536) * e2 ** 5 + (693693 / 1048576) * e2 ** 6 + (19324305 / 29360128) * e2 ** 7 + (4927697775 / 7516192768) * e2 ** 8
a2 = -(3 / 4) * e2 - (15 / 16) * e2 ** 2 - (525 / 512) * e2 ** 3 - (2205 / 2048) * e2 ** 4 - (72765 / 65536) * e2 ** 5 - (297297 / 262144) * e2 ** 6 - (135270135 / 117440512) * e2 ** 7 - (547521975 / 469762048) * e2 ** 8
a4 = (15 / 64) * e2 ** 2 + (105 / 256) * e2 ** 3 + (2205 / 4096) * e2 ** 4 + (10395 / 16384) * e2 ** 5 + (1486485 / 2097152) * e2 ** 6 + (45090045 / 58720256) * e2 **7 + (766530765 /939524096)*e**8
a6 = -(35/512)*e**3-(315/2048)*e**4-(31185/131072)*e**5-(165165/524288)*e**6-(45090045/117440512)*e**7-(209053845/335544320)*e**8
a8 = (315/16384)*e**4+(3465/65536)*e**5+(99099/1048576)*e**6+(4099095/29360128)*e**7+(348423075/1879048192)*e**8
a10 = -(693/131072)*e**5-(17325/262144)*e**6-(99099/524288)*e**7-(4099095/117440512)*e**8-(26801775/469762048)*e**9
a12 = (3003/2097152)*e**6+(315315/58720256)*e**7+(11486475/939524096)*e**8+(90810725/1879048192)*e**9
# 计算纬度的迭代初值
X1 = y/(a0*a)
Bf1 = X1+a1*math.sin(1*X1)+a3*math.sin(3*X1)+a5*math.sin(5*X1)+a7*math.sin(7*X1)+a9*math.sin(9*X1)+a11*math.sin(11*X1)
# 进行迭代计算，直到满足精度要求
while True:
Nf = a/math.sqrt(1-e*e*math.sin(Bf1)**)
tf = math.tan(Bf1)
Cf = ee*math.cos(Bf1)**)
Mf = a*((1-e*e/4-3*e*e*e*e/64-5*e*e*e*e*e*e/256)*Bf1-(3*e*e/8+3*e*e*e*e/32+45*e*e*e*e*e*e/1024)*math.sin(2*Bf1)+(15*e*e*e*e/256+45*e*e*e*e*e/1024)*math.sin(4*Bf1)-(35*e*e*e*e*e*e/3072)*math.sin(6*Bf1))
Bf = X1+(y-Mf)/Nf
if abs(Bf-Bf1) < 0.00000000001:
break
else:
Bf1 = Bf
# 计算经纬度坐标
lat_rad = Bf-tf*(x*x/(2*Nf*Nf)+x*x*x*x/(24*Nf*Nf*Nf*Nf)*(5+3*tf*tf+C*C-9*C*C*tf*tf)+x*x*x*x*x*x/(720*Nf*Nf*Nf*Nf*Nf*Nf)*(61+90*tf*tf+45*tf*tf*tf*tf))
lon_rad = Lon0Rad+x/(Nf*math.cos(Bf))-x*x*x/(6*Nf*Nf*Nf*math.cos(Bf))*(1+2*tf*tf+C*C)+x*x*x*x*x/(120*Nf*Nf*Nf*Nf*Nf*math.cos(Bf))*(5+28*tf*tf+24*tf*tf*tf*tf+6*C*C+8*C*C*tf*tf)
lon_deg = lon_rad / Rad
lat_deg = lat_rad / Rad
return lon_deg, lat_deg
# 定义一个函数，将BJS54经纬度坐标转换为WGS-84经纬度坐标
def BJS54ToWGS84(lon, lat):
# 将经纬度转换为弧度
lon_rad = lon * Rad
lat_rad = lat * Rad
# 将经纬度转换为XYZ坐标
N = a / math.sqrt(1 - e2 * math.sin(lat_rad) ** 2)
X = N * math.cos(lat_rad) * math.cos(lon_rad)
Y = N * math.cos(lat_rad) * math.sin(lon_rad)
Z = N * (1 - e2) * math.sin(lat_rad)
# 将XYZ坐标进行七参数转换，得到WGS-84下的XYZ坐标
X2 = (X - dx) / (1 + m) - rz * Y + ry * Z
Y2 = (Y - dy) / (1 + m) + rz * X - rx * Z
Z2 = (Z - dz) / (1 + m) - ry * X + rx * Y
# 将WGS-84下的XYZ坐标转换为经纬度坐标
lon2_rad = math.atan(Y2 / X2)
Bf = math.atan(Z2 / ((1 - e2) * math.sqrt(X2 ** 2 + Y2 ** 2)))
Nf = a / math.sqrt(1 - e2 * math.sin(Bf) ** 2)
lat2_rad = math.atan((Z2 + Nf * e2 * math.sin(Bf)) / math.sqrt(X2 ** 2 + Y2 ** 2))
lon2_deg = lon2_rad / Rad
lat2_deg = lat2_rad / Rad
return lon2_deg, lat2_deg
# 测试代码，将WGS-84下的北京市经纬度坐标转换为高斯-克吕格平面直角坐标
lon = 116.397128
lat = 39.916527
lon2, lat2 = WGS84ToBJS54(lon, lat)
x, y = BJS54ToGauss(lon2, lat2)
print("WGS-84经纬度坐标：", lon, lat)
print("BJS54经纬度坐标：", lon2, lat2)
print("高斯-克吕格平面直角坐标：", x, y)
# 输出结果
WGS-84经纬度坐标： 116.397128 39.916527
BJS54经纬度坐标： 116.39027611111111 39.913600000000004
高斯-克吕格平面直角坐标： 4428600.387 434801.886## 总结
本文介绍了地图坐标系的基本概念和原理，以及常见的地图坐标系的特点和应用场景。希望对有兴趣了解地图坐标系的读者有所帮助。如果您有任何疑问或建议，请在评论区留言。谢谢！
发表评论 取消回复