【精准三点定位求解汇总】利用Python或JavaScript高德地图开放平台实现精准三点定位（经纬度坐标与平面坐标转换法求解、几何绘图法求解）

众所周知，如果已知三个点的坐标，到一个未知点的距离，则可以利用以距离为半径画圆的方式来求得未知点坐标。
如果只有两个已知点，则只能得出两个未知点坐标，而第三个圆必定交于其中一个点

如图：
在这里插入图片描述
三个圆必定教于一个点

当然，如果第三个绿色圆圆心位于红蓝的圆心连线上，则依然交于两个点，所以在选择对照点时，应尽可能使对照点分布在未知点四周，多取几个点位未尝不是一门好事

几何绘图求解法（JavaScript）

之前发过一期这里再重新讲一遍

【优秀课设】JavaScript利用高德地图开放平台实现精准三点定位（通过距离半径，绘制多个circle）
blog.csdn.net/weixin_53403301/article/details/122329814

首先打开高德地图开发平台的覆盖物的添加与移除：
lbs.amap.com/demo/javascript-api/example/map-componets/map-overlays

其中，代码的第54行可以生成标记点
59-67可以画圆

注意两个主要参数：

center: new AMap.LngLat("116.403322", "39.920255"), // 圆心位置
radius: 1000,  //半径

半径r单位都是米

通过高德地图API：
lbs.amap.com/tools/picker
可以获取任意点位的坐标，只需要点击即可

比如武汉市坐标：114.304569,30.593354
在这里插入图片描述
将武汉市坐标复制到代码中，设置半径为300km，即300000

center: new AMap.LngLat("114.304569", "30.593354"), // 圆心位置
radius: 300000,  //半径

同时，定义circle2，添加第二个位置，如南昌市，设置半径为400km

var circle2 = new AMap.Circle({center: new AMap.LngLat("115.857972", "28.682976"), // 圆心位置radius: 400000,  //半径strokeColor: "#F33",  //线颜色strokeOpacity: 1,  //线透明度strokeWeight: 3,  //线粗细度fillColor: "#ee2200",  //填充颜色fillOpacity: 0.35 //填充透明度
});

再在后面加上语句，用于绘制circle和circle2

map.add(circle);
map.add(circle2);
map.setFitView();

点击右上角运行：
即可看到在地图上绘制了两个分别以武汉和南昌为半径的圆
在这里插入图片描述
两个圆的交点有两个，分别在常德附近和合肥附近
这里只用于演示，同样的方法可以添加第三个圆
如：

var circle3 = new AMap.Circle({center: new AMap.LngLat("116.378517", "39.865246"), // 圆心位置radius: 900000,  //半径strokeColor: "#F33",  //线颜色strokeOpacity: 1,  //线透明度strokeWeight: 3,  //线粗细度fillColor: "#ee2200",  //填充颜色fillOpacity: 0.35 //填充透明度
});
map.add(circle3);

在这里插入图片描述
59行开始，全部修改过的代码部分如下：

var circle = new AMap.Circle({center: new AMap.LngLat("114.304569", "30.593354"), // 圆心位置radius: 300000,  //半径strokeColor: "#F33",  //线颜色strokeOpacity: 1,  //线透明度strokeWeight: 3,  //线粗细度fillColor: "#ee2200",  //填充颜色fillOpacity: 0.35 //填充透明度
});var circle2 = new AMap.Circle({center: new AMap.LngLat("115.857972", "28.682976"), // 圆心位置radius: 400000,  //半径strokeColor: "#F33",  //线颜色strokeOpacity: 1,  //线透明度strokeWeight: 3,  //线粗细度fillColor: "#ee2200",  //填充颜色fillOpacity: 0.35 //填充透明度
});var circle3 = new AMap.Circle({center: new AMap.LngLat("116.378517", "39.865246"), // 圆心位置radius: 900000,  //半径strokeColor: "#F33",  //线颜色strokeOpacity: 1,  //线透明度strokeWeight: 3,  //线粗细度fillColor: "#ee2200",  //填充颜色fillOpacity: 0.35 //填充透明度
});
map.add(circle3);
map.add(circle);
map.add(circle2);
map.setFitView();

经纬度坐标与平面坐标转换求解法（Python）

主要用到这两个库：

import math
import pyproj

pyproj是用于坐标转换的这里采用的是utm平面坐标与WGS84经纬度坐标

这里的半径r单位都是米

高德地图用的是GCJ02坐标（还有腾讯） GPS用的是WGS84坐标（谷歌也是）百度地图用的是BD09坐标

所以在实际计算出来导入地图里面查看时要么采用WGS84地图要么就要涉及到坐标转换

首先是坐标转换
在utm坐标中有个zone值也就是经度分区计算公式为int(lon/6+31)
在反坐标转换中也要输入zone值这里直接可以输入转换前求得的zone值

def lonlat2utm(lon,lat):z=int(lon/6+31)proj = pyproj.Proj(proj='utm',zone=z,ellps='WGS84')return proj(lon, lat),zdef utm2lonlat(x,y,z):proj = pyproj.Proj(proj='utm',zone=z,ellps='WGS84')return proj(x, y,inverse=True)

平面上求圆的交点

def insec(p1,r1,p2,r2):x = p1[0]y = p1[1]R = r1a = p2[0]b = p2[1]S = r2d = math.sqrt((abs(a-x))**2 + (abs(b-y))**2)if d > (R+S) or d < (abs(R-S)):
#        print ("没有公共点")return elif d == 0 and R==S :
#        print ("两个圆同心")returnelse:A = (R**2 - S**2 + d**2) / (2 * d)h = math.sqrt(R**2 - A**2)x2 = x + A * (a-x)/dy2 = y + A * (b-y)/dx3 = round(x2 - h * (b - y) / d,2)y3 = round(y2 + h * (a - x) / d,2)x4 = round(x2 + h * (b - y) / d,2)y4 = round(y2 - h * (a - x) / d,2)c1=[x3, y3]c2=[x4, y4]return c1,c2

经纬度上求两个圆交点坐标
在反坐标转换中也要输入zone值这里直接可以输入转换前求得的两个经纬度坐标的zone平均值
所以两个要求的区域离得不远误差就很小离得太远了误差就可能大需要手动去调整这里是通用函数

def location_trans(p1,r1,p2,r2):z1=lonlat2utm(p1[0],p1[1])z2=lonlat2utm(p2[0],p2[1])z=int((z1[1]+z2[1])/2)    C=insec(z1[0],r1,z2[0],r2)if C:a=utm2lonlat(C[0][0],C[0][1],z)b=utm2lonlat(C[1][0],C[1][1],z)c1=[a[0], a[1]]c2=[b[0], b[1]]return c1,c2else:return

运行：

 a=[[114.304569,30.593354],300000]b=[[115.857972,28.682976],400000]c=[[116.378517,39.865246],900000]print(location_trans(b[0],b[1],c[0],c[1]))

输出：

([114.12482189881902, 31.962881802790577], [117.87031764680636, 31.841927527011755])

在求三个圆的交点时最多会求出六个点
在六个点中筛选出离另外一个圆最近的点即可得出三个相近点的坐标

求离得近的那个点的平面坐标

def location_min(p1,p2,p,r):d1=math.fabs(r-math.sqrt((p[0]-p1[0])**2+(p[1]-p1[1])**2))d2=math.fabs(r-math.sqrt((p[0]-p2[0])**2+(p[1]-p2[1])**2))if d1<d2:return p1else:return p2

得到三个点后求三个点的中心点即可算出大概位置
中心点的zone值是根据其他三个zone值求平均来确定的
所以三个要求的区域离得不远误差就很小离得太远了误差就可能大需要手动去调整这里是通用函数

def location_judg(p1,r1,p2,r2,p3,r3):li=[]z1=lonlat2utm(p1[0],p1[1])z2=lonlat2utm(p2[0],p2[1])z3=lonlat2utm(p3[0],p3[1])z12=int((z1[1]+z2[1])/2)z13=int((z1[1]+z3[1])/2)z23=int((z2[1]+z3[1])/2)z=int((z12+z13+z23)/3)C12=insec(z1[0],r1,z2[0],r2)C13=insec(z1[0],r1,z3[0],r3)C23=insec(z2[0],r2,z3[0],r3)if C12:m12=location_min(C12[0],C12[1],z3[0],r3)li.append(utm2lonlat(m12[0],m12[1],z12))else:li.append(None)if C13:m13=location_min(C13[0],C13[1],z2[0],r2)li.append(utm2lonlat(m13[0],m13[1],z13))else:li.append(None)if C23:m23=location_min(C23[0],C23[1],z1[0],r1)li.append(utm2lonlat(m23[0],m23[1],z23))else:li.append(None)if C12 and C13 and C23:
#        print("三个坐标作的圆都有公共点")m=[(m12[0]+m13[0]+m23[0])/3,(m12[1]+m13[1]+m23[1])/3]li.append(utm2lonlat(m[0],m[1],z))return lielif C12 or C13 or C23:
#        print("三个坐标作的圆不全有公共点")li.append(None)return lielse:
#        print("三个坐标作的圆都没有公共点")return

最后返回的列表分别是12的最接近坐标 13的最接近坐标 23的最接近坐标和这三个坐标的中心点坐标如果不存在则返回None

运行：

if __name__ == "__main__":a=[[114.304569,30.593354],300000]b=[[115.857972,28.682976],400000]c=[[116.378517,39.865246],900000]print(location_trans(b[0],b[1],c[0],c[1]))    print(location_judg(a[0],a[1],b[0],b[1],c[0],c[1]))

结果：

[(116.85351953263574, 32.18782636821823), (117.13697531307241, 31.774218803048125), (117.87031764680636, 31.841927527011755), (117.28744847106574, 31.935380071325877)]

整体代码：

# -*- coding: utf-8 -*-
import math
import pyprojdef lonlat2utm(lon,lat):z=int(lon/6+31)proj = pyproj.Proj(proj='utm',zone=z,ellps='WGS84')return proj(lon, lat),zdef utm2lonlat(x,y,z):proj = pyproj.Proj(proj='utm',zone=z,ellps='WGS84')return proj(x, y,inverse=True)def insec(p1,r1,p2,r2):x = p1[0]y = p1[1]R = r1a = p2[0]b = p2[1]S = r2d = math.sqrt((abs(a-x))**2 + (abs(b-y))**2)if d > (R+S) or d < (abs(R-S)):
#        print ("没有公共点")return elif d == 0 and R==S :
#        print ("两个圆同心")returnelse:A = (R**2 - S**2 + d**2) / (2 * d)h = math.sqrt(R**2 - A**2)x2 = x + A * (a-x)/dy2 = y + A * (b-y)/dx3 = round(x2 - h * (b - y) / d,2)y3 = round(y2 + h * (a - x) / d,2)x4 = round(x2 + h * (b - y) / d,2)y4 = round(y2 - h * (a - x) / d,2)c1=[x3, y3]c2=[x4, y4]return c1,c2def location_trans(p1,r1,p2,r2):z1=lonlat2utm(p1[0],p1[1])z2=lonlat2utm(p2[0],p2[1])z=int((z1[1]+z2[1])/2)    C=insec(z1[0],r1,z2[0],r2)if C:a=utm2lonlat(C[0][0],C[0][1],z)b=utm2lonlat(C[1][0],C[1][1],z)c1=[a[0], a[1]]c2=[b[0], b[1]]return c1,c2else:returndef location_min(p1,p2,p,r):d1=math.fabs(r-math.sqrt((p[0]-p1[0])**2+(p[1]-p1[1])**2))d2=math.fabs(r-math.sqrt((p[0]-p2[0])**2+(p[1]-p2[1])**2))if d1<d2:return p1else:return p2def location_judg(p1,r1,p2,r2,p3,r3):li=[]z1=lonlat2utm(p1[0],p1[1])z2=lonlat2utm(p2[0],p2[1])z3=lonlat2utm(p3[0],p3[1])z12=int((z1[1]+z2[1])/2)z13=int((z1[1]+z3[1])/2)z23=int((z2[1]+z3[1])/2)z=int((z12+z13+z23)/3)C12=insec(z1[0],r1,z2[0],r2)C13=insec(z1[0],r1,z3[0],r3)C23=insec(z2[0],r2,z3[0],r3)if C12:m12=location_min(C12[0],C12[1],z3[0],r3)li.append(utm2lonlat(m12[0],m12[1],z12))else:li.append(None)if C13:m13=location_min(C13[0],C13[1],z2[0],r2)li.append(utm2lonlat(m13[0],m13[1],z13))else:li.append(None)if C23:m23=location_min(C23[0],C23[1],z1[0],r1)li.append(utm2lonlat(m23[0],m23[1],z23))else:li.append(None)if C12 and C13 and C23:
#        print("三个坐标作的圆都有公共点")m=[(m12[0]+m13[0]+m23[0])/3,(m12[1]+m13[1]+m23[1])/3]li.append(utm2lonlat(m[0],m[1],z))return lielif C12 or C13 or C23:
#        print("三个坐标作的圆不全有公共点")li.append(None)return lielse:
#        print("三个坐标作的圆都没有公共点")return if __name__ == "__main__":a=[[114.304569,30.593354],300000]b=[[115.857972,28.682976],400000]c=[[116.378517,39.865246],900000]print(location_trans(b[0],b[1],c[0],c[1]))    print(location_judg(a[0],a[1],b[0],b[1],c[0],c[1]))