我想找出一个未知的目标位置(纬度和经度坐标)。有3个已知点(纬度和经度坐标对),每个点到目标位置的距离(以千米为单位)。如何计算目标位置的坐标?

例如,说我有以下数据点

37.418436,-121.963477   0.265710701754km
37.417243,-121.961889   0.234592423446km
37.418692,-121.960194   0.0548954278262km


我要做什么我想知道一个函数的数学运算,该函数将其作为输入并返回37.417959,-121.961954作为输出。 .movable-type.co.uk / scripts / latlong.html我了解这样的基本原理:有了这样的三个圆圈,您恰好得到一个重叠点。我不知道的是使用此输入来计算该点所需的数学。

评论

这是一个页面,引导您完成寻找三个坐标中心的数学过程。也许这可能会有所帮助。

这是否需要在球体/椭球面上,还是可以使用平面算法?

我无法给您答案,但我想我可以为您指明正确的方向。三个坐标=三个中心点。三个距离=三个圆圈。相交的两个圆可以有一个/一个/两个解的可能性。三个圆可以没有/一个/或一个区域作为其解决方案。获得三个圆的圆公式,并用方程组/代数系统求解。

实际上,您甚至不需要系统来解决这一问题。有一种或两种可能性,但是由于您具有距离值,因此可以分离出正确的答案。

+1这是一个好问题。起初,我认为可以通过google轻松找到解决方案,但显然不是。也许可以更笼统地说出问题:给定N个点,每个点不仅具有距离,而且具有误差范围,找到置信椭圆。

#1 楼

在浏览了一下Wikipedia并在StackOverflow上找到了相同的问题/答案之后,我想我会对此进行尝试,并尝试填补空白。输出,但似乎是错误的。我在ArcMap中绘制了点,将其缓冲到指定的距离,在缓冲区上相交到,然后捕获相交的顶点以获得解。您建议的输出是绿色的点。我计算了标注框中的值,该值大约是ArcMap给出的用于从相交得出的解决方案的值的3米。太糟糕了,只需要将您的大地坐标转换为笛卡尔ECEF,即可在此处找到。如果您不使用椭圆体,则a / x + h项可以替换为真实球体半径。在python

 import math
import numpy

#assuming elevation = 0
earthR = 6371
LatA = 37.418436
LonA = -121.963477
DistA = 0.265710701754
LatB = 37.417243
LonB = -121.961889
DistB = 0.234592423446
LatC = 37.418692
LonC = -121.960194
DistC = 0.0548954278262

#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
#   1. Convert Lat/Long to radians
#   2. Convert Lat/Long(radians) to ECEF
xA = earthR *(math.cos(math.radians(LatA)) * math.cos(math.radians(LonA)))
yA = earthR *(math.cos(math.radians(LatA)) * math.sin(math.radians(LonA)))
zA = earthR *(math.sin(math.radians(LatA)))

xB = earthR *(math.cos(math.radians(LatB)) * math.cos(math.radians(LonB)))
yB = earthR *(math.cos(math.radians(LatB)) * math.sin(math.radians(LonB)))
zB = earthR *(math.sin(math.radians(LatB)))

xC = earthR *(math.cos(math.radians(LatC)) * math.cos(math.radians(LonC)))
yC = earthR *(math.cos(math.radians(LatC)) * math.sin(math.radians(LonC)))
zC = earthR *(math.sin(math.radians(LatC)))

P1 = numpy.array([xA, yA, zA])
P2 = numpy.array([xB, yB, zB])
P3 = numpy.array([xC, yC, zC])

#from wikipedia
#transform to get circle 1 at origin
#transform to get circle 2 on x axis
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = numpy.dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = numpy.dot(ey, P3 - P1)

#from wikipedia
#plug and chug using above values
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x)

# only one case shown here
z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))

#triPt is an array with ECEF x,y,z of trilateration point
triPt = P1 + x*ex + y*ey + z*ez

#convert back to lat/long from ECEF
#convert to degrees
lat = math.degrees(math.asin(triPt[2] / earthR))
lon = math.degrees(math.atan2(triPt[1],triPt[0]))

print lat, lon
 


评论


我本来打算给出一个类似的答案,但是现在没有必要了!得到我的支持。

–拉斯
2010年7月28日在15:54

麻木的救援!当将“ triPt”替换为“ triLatPt”时,它将编译,但否则返回37.4191023738 -121.960579208。做得好

–WolfOdrade
2010年7月28日在20:18

干得好!如果我将地理坐标系替换为本地[Cartesian]坐标系,这是否仍然有效?

– Zengr
13年10月3日,19:47

对于那些在c ++域中的人。

–raaj
2014年8月14日17:32



谢谢@wwnick!我已将其移植到JavaScript(用于Node,但可以轻松转换为在浏览器中工作)。 gist.github.com/dav-/bb7103008cdf9359887f

– HellaMad
2014年11月26日,下午1:34

#2 楼

我不确定我是否还很幼稚,但是,如果按大小缓冲每个点,然后将所有三个圆相交会得到正确的位置?

可以使用以下方法计算交点:空间API。示例:


GeoScript
Java拓扑套件
NET拓扑套件
GEOS


评论


确实,他对获取该交点的公式感兴趣。

– Vinko Vrsalovic
2010年7月22日在20:51

使用空间API,无需使用纯数学就可以做到这一点。

–乔治·席尔瓦(George Silva)
2010年7月23日下午0:31

@George您能举一个这样的API的例子吗?

– nohat
2010年7月23日下午2:00

编辑帖子以反映nohat的要求。

–乔治·席尔瓦(George Silva)
2010年7月23日在2:23

+1,具有良好的横向思维能力,即使计算效率最高!

–fmark
2010年7月23日在8:17

#3 楼

以下注释使用了等平面几何(即,您必须将坐标投影到适当的局部坐标系中)。

我的推理,下面是使用Python的实例: >获取2个数据点(分别称为ab)。称为目标点x。我们已经知道距离axbx。我们可以使用毕达哥拉斯定理计算距离ab

>>> import math
>>> a = (1, 4)
>>> b = (3, 6)
>>> dist_ax = 3
>>> dist_bx = 5.385
# Pythagoras's theorem
>>> dist_ab = math.sqrt(abs(a[0]-b[0])**2 + abs(a[1]-b[1])**2)
>>> dist_ab
2.8284271247461903


现在,您可以算出这些线的角度:

>>> angle_abx = math.acos((dist_bx * dist_bx + dist_ab * dist_ab - dist_ax * dist_ax)/(2 * dist_bx * dist_ab))
>>> math.degrees(angle_abx)
23.202973815040256
>>> angle_bax = math.acos((dist_ax * dist_ax + dist_ab * dist_ab - dist_bx * dist_bx)/(2 * dist_ax * dist_ab))
>>> math.degrees(angle_bax)
134.9915256259537
>>> angle_axb = math.acos((dist_ax * dist_ax + dist_bx * dist_bx - dist_ab * dist_ab)/(2 * dist_ax * dist_bx))
>>> math.degrees(angle_axb)
21.805500559006095


很不幸,我来不及为您完成答案,但是,现在您知道了角度,您可以计算x的两个可能位置。然后,使用第三点c可以计算出正确的位置。

#4 楼

这可能有效。在python中再次快速地,您可以将其放在函数主体中xN,yN =点的坐标,r1和r2 =半径值

dX = x2 - x1
dY = y2 - y1

centroidDistance = math.sqrt(math.pow(e,2) + math.pow(dY,2)) #distance from centroids
distancePL = (math.pow(centroidDistance,2) + (math.pow(r1,2) - math.pow(r2,2))) / (2 * centroidDistance) #distance from point to a line splitting the two centroids

rx1 = x1 + (dX *k)/centroidDistance + (dY/centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))
ry1 = y1 + (dY*k)/centroidDistance - (dX /centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))

rx2 = x1 + (dX *k)/centroidDistance - (dY/centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))
ry2 = y1 + (dY*k)/centroidDistance + (dX /centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))


rx和ry值是一个圆上两个交点的返回值(应该在数组中),如果这有助于弄清楚事情的话。如果第一次迭代的结果与第二次迭代的结果相比较(可能在一定公差范围内),那么您就有了交点。这不是一个很好的解决方案,尤其是当您开始在过程中添加多个点时,但这是最简单的方法,无需解决方程组。

评论


您的代码中的“ e”和“ k”是什么?

– ReinierDG
2013年9月14日上午8:09

我不记得了:-) wwnick的答案更多的是如果您只有三个圆圈,您想实现的目标。

–WolfOdrade
2013年9月16日15:15

#5 楼

您可以使用postgis中的空间API(St_Intersection,St_buffer函数)。
正如fmark所注意到的,您还必须记住Postgis使用平面算法,但是对于小区域,使用等距投影不会带来太多错误。 />

评论


PostGIS可以使用GEOGRAPHY类型而不是GEOMETRY类型进行球体计算。

–fmark
2010年7月25日在22:47

#6 楼

用PHP语言执行:

//assuming elevation = 0
$earthR = 6371; // in km ( = 3959 in miles)

$LatA = 37.418436;
$LonA = -121.963477;
$DistA = 0.265710701754;

$LatB = 37.417243;
$LonB = -121.961889;
$DistB = 0.234592423446;

$LatC = 37.418692;
$LonC = -121.960194;
$DistC = 0.0548954278262;

/*
#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
#   1. Convert Lat/Long to radians
#   2. Convert Lat/Long(radians) to ECEF
*/
$xA = $earthR *(cos(deg2rad($LatA)) * cos(deg2rad($LonA)));
$yA = $earthR *(cos(deg2rad($LatA)) * sin(deg2rad($LonA)));
$zA = $earthR *(sin(deg2rad($LatA)));

$xB = $earthR *(cos(deg2rad($LatB)) * cos(deg2rad($LonB)));
$yB = $earthR *(cos(deg2rad($LatB)) * sin(deg2rad($LonB)));
$zB = $earthR *(sin(deg2rad($LatB)));

$xC = $earthR *(cos(deg2rad($LatC)) * cos(deg2rad($LonC)));
$yC = $earthR *(cos(deg2rad($LatC)) * sin(deg2rad($LonC)));
$zC = $earthR *(sin(deg2rad($LatC)));

/*
INSTALL:
sudo pear install Math_Vector-0.7.0
sudo pear install Math_Matrix-0.8.7
*/
// Include PEAR::Math_Matrix
//  /usr/share/php/Math/Matrix.php
//  include_path=".:/usr/local/php/pear/"
require_once 'Math/Matrix.php';
require_once 'Math/Vector.php';
require_once 'Math/Vector3.php';


$P1vector = new Math_Vector3(array($xA,$yA,$zA));
$P2vector = new Math_Vector3(array($xB,$yB,$zB));
$P3vector = new Math_Vector3(array($xC,$yC,$zC));

#from wikipedia: http://en.wikipedia.org/wiki/Trilateration
#transform to get circle 1 at origin
#transform to get circle 2 on x axis

// CALC EX
$P2minusP1 = Math_VectorOp::substract($P2vector, $P1vector);
$l = new Math_Vector($P2minusP1);
$P2minusP1_length = $l->length();
$norm = new Math_Vector3(array($P2minusP1_length,$P2minusP1_length,$P2minusP1_length));
$d = $norm; //save calc D
$ex = Math_VectorOp::divide($P2minusP1, $norm);
//echo "ex: ".$ex->toString()."\n";
$ex_x = floatval( $ex->_tuple->getData()[0] ) ;
$ex_y = floatval( $ex->_tuple->getData()[1] ) ;
$ex_z = floatval( $ex->_tuple->getData()[2] ) ;
$ex = new Math_Vector3(array($ex_x,$ex_y,$ex_z));

// CALC i
$P3minusP1 = Math_VectorOp::substract($P3vector, $P1vector);
$P3minusP1_x = floatval( $P3minusP1->_tuple->getData()[0] ) ;
$P3minusP1_y = floatval( $P3minusP1->_tuple->getData()[1] ) ;
$P3minusP1_z = floatval( $P3minusP1->_tuple->getData()[2] ) ;
$P3minusP1 = new Math_Vector3(array($P3minusP1_x,$P3minusP1_y,$P3minusP1_z));
$i = Math_VectorOp::dotProduct($ex, $P3minusP1);
//echo "i = $i\n";

// CALC EY
$iex = Math_VectorOp::scale($i, $ex);
//echo " iex = ".$iex->toString()."\n";
$P3P1iex = Math_VectorOp::substract($P3minusP1, $iex);
//echo " P3P1iex = ".$P3P1iex->toString()."\n";
$l = new Math_Vector($P3P1iex);
$P3P1iex_length = $l->length();
$norm = new Math_Vector3(array($P3P1iex_length,$P3P1iex_length,$P3P1iex_length));
//echo "norm: ".$norm->toString()."\n";
$ey = Math_VectorOp::divide($P3P1iex, $norm);
//echo " ey = ".$ey->toString()."\n";
$ey_x = floatval( $ey->_tuple->getData()[0] ) ;
$ey_y = floatval( $ey->_tuple->getData()[1] ) ;
$ey_z = floatval( $ey->_tuple->getData()[2] ) ;
$ey = new Math_Vector3(array($ey_x,$ey_y,$ey_z));

// CALC EZ
$ez = Math_VectorOp::crossProduct($ex, $ey);
//echo " ez = ".$ez->toString()."\n";

// CALC D
// do it before
$d = floatval( $d->_tuple->getData()[0] ) ;
//echo "d = $d\n";

// CALC J
$j = Math_VectorOp::dotProduct($ey, $P3minusP1);
//echo "j = $j\n";

#from wikipedia
#plug and chug using above values
$x = (pow($DistA,2) - pow($DistB,2) + pow($d,2))/(2*$d);
$y = ((pow($DistA,2) - pow($DistC,2) + pow($i,2) + pow($j,2))/(2*$j)) - (($i/$j)*$x);

# only one case shown here
$z = sqrt( pow($DistA,2) - pow($x,2) - pow($y,2) );

//echo "x = $x - y = $y  - z = $z\n";

#triPt is an array with ECEF x,y,z of trilateration point
$xex = Math_VectorOp::scale($x, $ex);
$yey = Math_VectorOp::scale($y, $ey);
$zez = Math_VectorOp::scale($z, $ez);

// CALC $triPt = $P1vector + $xex + $yey + $zez;
$triPt = Math_VectorOp::add($P1vector, $xex);
$triPt = Math_VectorOp::add($triPt, $yey);
$triPt = Math_VectorOp::add($triPt, $zez);
//echo " triPt = ".$triPt->toString()."\n";
$triPt_x = floatval( $triPt->_tuple->getData()[0] ) ;
$triPt_y = floatval( $triPt->_tuple->getData()[1] ) ;
$triPt_z = floatval( $triPt->_tuple->getData()[2] ) ;


#convert back to lat/long from ECEF
#convert to degrees
$lat = rad2deg(asin($triPt_z / $earthR));
$lon = rad2deg(atan2($triPt_y,$triPt_x));

echo $lat .','. $lon;