这似乎与这个问题相反(纬度/经度点之间的距离)。
我已经研究过Haversine公式,并认为它与世界的近似值已经足够接近了。
我假设我需要解决未知数的Haversine公式经/纬度,这是正确的吗?有没有很好的网站谈论这种事情?看来这很常见,但是我的谷歌搜索只出现了与上述问题类似的问题。
我真正想要的只是一个公式。我想给它起一个纬度/经度,方位角和距离(英里或公里),我想从中得出一个纬度/经度对,代表如果他们沿途旅行,最终会去哪里那条路线。
#1 楼
我很好奇这个公式的结果与Esri的pe.dll相比如何。(引用)。
点{lat,lon}是一个距点1的tc半径上的距离d out
,如果:
lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
IF (cos(lat)=0)
lon=lon1 // endpoint a pole
ELSE
lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
ENDIF
该算法限于距离
使得dlon
经度上延伸。如果允许更大的距离,则需要一个完全通用但更复杂的算法
:
lat =asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
dlon=atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(lat))
lon=mod( lon1-dlon +pi,2*pi )-pi
这里是用于测试的html页面。
评论
感谢你能这么快回复。让我消化一些信息,我会与您联系。从表面上看,这看起来确实不错。
–Jason Whitehorn
2011-2-4在20:21
在从html页面检索32N,117W到40N,82W的距离和正向方位角后,我尝试使用pe.dll(实际上是solaris上的libpe.so)尝试了直接案例。使用32N,117W,d = 3255.056515890041,azi = 64.24498012065699,我得到40N,82W(实际上是82.00000000064)。
– mkennedy
2011年2月4日在22:17
太棒了!非常感谢您与Ed Williams撰写的《航空配方》文章的链接,我之前从未见过,但到目前为止,它已被证明是一本好书。只是对将来查看此内容的任何人的注释,此公式的输入和输出都是弧度,甚至是距离。
–Jason Whitehorn
2011-2-5在6:09
此公式中的距离单位是多少?
– Karthik Murugan
2015年4月2日在6:35
@KarthikMurugan Ed的介绍说,距离单位沿大圆弧度是弧度。
– Kirk Kuykendall
17年9月6日在18:35
#2 楼
如果您在飞机上,那么在北向东度数方位处r米处的点在北向会错开r * cos(a),在东向会错开r * sin(a)。 (这些语句或多或少定义了正弦和余弦。)尽管您不在平面中,但您正在模拟地球表面的弯曲椭圆体的表面上任意距离不到几百公里的地方覆盖了表面的一小部分,因此在大多数实际应用中可以认为是平坦的。唯一剩下的复杂性是,一个经度不能覆盖与一个纬度相同的距离。在球形地球模型中,一个经度只有cos(纬度)与纬度一样长。 (在椭圆模型中,这仍然是一个很好的近似值,大约为2.5个有效数字。)
最后,一个纬度约为10 ^ 7/90 = 111,111米。现在,我们拥有将米转换为度所需的所有信息:
向北位移为r * cos(a)/ 111111度;
向东位移为r * sin (a)/ cos(纬度)/ 111111度。
例如,在-0.31399度的纬度和北纬30度的方位角上,我们可以计算
cos(a) = cos(30 degrees) = cos(pi/6 radians) = Sqrt(3)/2 = 0.866025.
sin(a) = sin(30 degrees) = sin(pi/6 radians) = 1/2 = 0.5.
cos(latitude) = cos(-0.31399 degrees) = cos(-0.00548016 radian) = 0.999984984.
r = 100 meters.
east displacement = 100 * 0.5 / 0.999984984 / 111111 = 0.000450007 degree.
north displacement = 100 * 0.866025 / 111111 = 0.000779423 degree.
从(-78.4437,-0.31399)开始,新位置在(-78.4437 + 0.00045,-0.31399 + 0.0007794)=(-78.4432,-0.313211)。 br />
在现代ITRF00参考系统中,一个更准确的答案是(-78.4433,-0.313207):与近似答案相距0.43米,表明在这种情况下的近似误差为0.43%。要获得更高的精度,您必须使用椭圆距离公式(复杂得多)或零散度的高保真共形投影(以便使轴承正确)。
评论
+1以正确理解数学上下文(即其局部平面)
–弗兰克·康利(Frank Conry)
2014年8月10日下午5:46
#3 楼
如果您需要JavaScript解决方案,请考虑以下这些functions
和这个小提琴:var gis = {
/**
* All coordinates expected EPSG:4326
* @param {Array} start Expected [lon, lat]
* @param {Array} end Expected [lon, lat]
* @return {number} Distance - meter.
*/
calculateDistance: function(start, end) {
var lat1 = parseFloat(start[1]),
lon1 = parseFloat(start[0]),
lat2 = parseFloat(end[1]),
lon2 = parseFloat(end[0]);
return gis.sphericalCosinus(lat1, lon1, lat2, lon2);
},
/**
* All coordinates expected EPSG:4326
* @param {number} lat1 Start Latitude
* @param {number} lon1 Start Longitude
* @param {number} lat2 End Latitude
* @param {number} lon2 End Longitude
* @return {number} Distance - meters.
*/
sphericalCosinus: function(lat1, lon1, lat2, lon2) {
var radius = 6371e3; // meters
var dLon = gis.toRad(lon2 - lon1),
lat1 = gis.toRad(lat1),
lat2 = gis.toRad(lat2),
distance = Math.acos(Math.sin(lat1) * Math.sin(lat2) +
Math.cos(lat1) * Math.cos(lat2) * Math.cos(dLon)) * radius;
return distance;
},
/**
* @param {Array} coord Expected [lon, lat] EPSG:4326
* @param {number} bearing Bearing in degrees
* @param {number} distance Distance in meters
* @return {Array} Lon-lat coordinate.
*/
createCoord: function(coord, bearing, distance) {
/** http://www.movable-type.co.uk/scripts/latlong.html
* φ is latitude, λ is longitude,
* θ is the bearing (clockwise from north),
* δ is the angular distance d/R;
* d being the distance travelled, R the earth’s radius*
**/
var
radius = 6371e3, // meters
δ = Number(distance) / radius, // angular distance in radians
θ = gis.toRad(Number(bearing));
φ1 = gis.toRad(coord[1]),
λ1 = gis.toRad(coord[0]);
var φ2 = Math.asin(Math.sin(φ1)*Math.cos(δ) +
Math.cos(φ1)*Math.sin(δ)*Math.cos(θ));
var λ2 = λ1 + Math.atan2(Math.sin(θ) * Math.sin(δ)*Math.cos(φ1),
Math.cos(δ)-Math.sin(φ1)*Math.sin(φ2));
// normalise to -180..+180°
λ2 = (λ2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI;
return [gis.toDeg(λ2), gis.toDeg(φ2)];
},
/**
* All coordinates expected EPSG:4326
* @param {Array} start Expected [lon, lat]
* @param {Array} end Expected [lon, lat]
* @return {number} Bearing in degrees.
*/
getBearing: function(start, end){
var
startLat = gis.toRad(start[1]),
startLong = gis.toRad(start[0]),
endLat = gis.toRad(end[1]),
endLong = gis.toRad(end[0]),
dLong = endLong - startLong;
var dPhi = Math.log(Math.tan(endLat/2.0 + Math.PI/4.0) /
Math.tan(startLat/2.0 + Math.PI/4.0));
if (Math.abs(dLong) > Math.PI) {
dLong = (dLong > 0.0) ? -(2.0 * Math.PI - dLong) : (2.0 * Math.PI + dLong);
}
return (gis.toDeg(Math.atan2(dLong, dPhi)) + 360.0) % 360.0;
},
toDeg: function(n) { return n * 180 / Math.PI; },
toRad: function(n) { return n * Math.PI / 180; }
};
因此,如果要计算新坐标,可以这样:
var start = [15, 38.70250];
var end = [21.54967, 38.70250];
var total_distance = gis.calculateDistance(start, end); // meters
var percent = 10;
// this can be also meters
var distance = (percent / 100) * total_distance;
var bearing = gis.getBearing(start, end);
var new_coord = gis.createCoord(icon_coord, bearing, distance);
评论
在php中也是如此:gist.github.com/AidasK/17311a6a32950954335144e8977bf8cd
–艾达斯
3月1日13:48
#4 楼
我在ObjectiveC中完成了这项工作。这里的关键是要知道您需要以弧度为单位的经/纬度点,并且需要在应用方程式之后将它们转换回经/纬度。另外,要知道距离和tc都以弧度为单位。这是原始公式:
lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
IF (cos(lat)=0)
lon=lon1 // endpoint a pole
ELSE
lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
ENDIF
在ObjC中实现,其中弧度是从N逆时针测量的弧度(例如PI / 2为W,PI为S,2 PI / 3为E),距离以公里为单位。
+ (CLLocationCoordinate2D)displaceLatLng:(CLLocationCoordinate2D)coordinate2D withRadian:(double)radian
withDistance:(CGFloat)distance {
double lat1Radians = coordinate2D.latitude * (M_PI / 180);
double lon1Radians = coordinate2D.longitude * (M_PI / 180);
double distanceRadians = distance / 6371;
double lat = asin(sin(lat1Radians)*cos(distanceRadians)+cos(lat1Radians)*sin(distanceRadians)
*cos(radian));
double lon;
if (cos(lat) == 0) {
lon = lon1Radians;
} else {
lon = fmodf((float) (lon1Radians - asin(sin(radian) * sin(distanceRadians) / cos(lat)) + M_PI),
(float) (2 * M_PI)) - M_PI;
}
double newLat = lat * (180 / M_PI);
double newLon = lon * (180 / M_PI);
return CLLocationCoordinate2DMake(newLat, newLon);
}
评论
我正在寻找一种解决方案,希望从我站立的位置到北50英里,西50英里,东50英里等等的地方得到4纬度,长...在上面的答案中,您可以解释我如何实现这个 ?
–拉胡尔·维斯(Rahul Vyas)
18-2-20在12:29
#5 楼
如果您对提高准确性感兴趣,可以找Vincenty。 (链接指向的是“直接”形式,这正是您所追求的。)如果有代码,则有很多现有的实现。
还有一个问题:您不是在假设旅行者在整个旅程中都保持相同的方位,对吗?如果是这样,那么这些方法将无法回答正确的问题。 (最好投影到墨卡托,画一条直线,然后取消投影结果。)
评论
很好的问题,尽管我的问题中的措辞可能表明我正在为旅行者计算目的地,但事实并非如此。好点。这样做主要是为了让我可以计算出边界区域(小范围,例如50英里)……我希望这是有道理的。
–Jason Whitehorn
2011年2月6日,下午3:45
gis.stackexchange.com/questions/3264/…也有同样的问题(从点和距离构造边界区域)
– Dan S.
2011-2-7在19:46
#6 楼
这是Python解决方案: def displace(self, theta, distance):
"""
Displace a LatLng theta degrees counterclockwise and some
meters in that direction.
Notes:
http://www.movable-type.co.uk/scripts/latlong.html
0 DEGREES IS THE VERTICAL Y AXIS! IMPORTANT!
Args:
theta: A number in degrees.
distance: A number in meters.
Returns:
A new LatLng.
"""
theta = np.float32(theta)
delta = np.divide(np.float32(distance), np.float32(E_RADIUS))
def to_radians(theta):
return np.divide(np.dot(theta, np.pi), np.float32(180.0))
def to_degrees(theta):
return np.divide(np.dot(theta, np.float32(180.0)), np.pi)
theta = to_radians(theta)
lat1 = to_radians(self.lat)
lng1 = to_radians(self.lng)
lat2 = np.arcsin( np.sin(lat1) * np.cos(delta) +
np.cos(lat1) * np.sin(delta) * np.cos(theta) )
lng2 = lng1 + np.arctan2( np.sin(theta) * np.sin(delta) * np.cos(lat1),
np.cos(delta) - np.sin(lat1) * np.sin(lat2))
lng2 = (lng2 + 3 * np.pi) % (2 * np.pi) - np.pi
return LatLng(to_degrees(lat2), to_degrees(lng2))
#7 楼
在给定方位角和与上一个坐标的距离的情况下,我使用下面描述的方法确定下一个坐标。我从互联网上读取的其他方法在准确性方面存在问题。我用它来确定土地的面积,它是一个多边形,并在Google Earth中绘制该多边形。土地标题以这种方式写有方位角和距离:“北或南x度y分钟东或西,z米到n点;”。
因此,从参考点的给定坐标开始,我首先使用Haversine公式计算每1度纬度和1度经度的距离,因为这会根据位置而变化。然后,我从三角正弦和余弦公式确定下一个坐标。
下面是javascript:
var mapCenter = new google.maps.LatLng(referencePointLatitude, referencePointLongitude); //the ref point lat and lon must be given, usually a land mark (BLLM)
var latDiv = latDiv(mapCenter); //distance per one degree latitude in this location
var lngDiv = lngDiv(mapCenter); //distance per one degree longitude in this location
var LatLng2 = NextCoord(PrevCoord,NorthOrSouth,x,y,EastOrWest,z); //next coordinate given the bearing and distance from previous coordinate
var Lat2 = LatLng2.lat(); //next coord latitude in degrees
var Lng2 = LatLng2.lng(); //next coord longitude in degrees
var polygon=[p1,p2,...,pn-1,pn,p1]; //p1,p2,etc. are the coordinates of points of the polygon, i.e. the land Title. Be sure to close the polygon to the point of beginning p1
var area = Area(polygon); //area of the polygon in sq.m.
function NextCoord(PrevCoord,NorthOrSouth,x,y,EastOrWest,z) {
var angle = ( x + ( y / 60 ) ) * Math.PI / 180;
var a = 1;
var b = 1;
if (NorthOrSouth == 'South') { a = -1; }
if (EastOrWest == 'West') { b = -1; }
var nextLat = PrevCoord.lat() + ( a * z * Math.cos(angle) / latDiv );
var nextLng = PrevCoord.lng() + ( b * z * Math.sin(angle) / lngDiv );
var nextCoord = new google.maps.LatLng(nextLat, nextLng);
return nextCoord;
}
function latDiv(mapCenter) {
var p1 = new google.maps.LatLng(mapCenter.lat()-0.5, mapCenter.lng());
var p2 = new google.maps.LatLng(mapCenter.lat()+0.5, mapCenter.lng());
return dist(p1,p2);
}
function lngDiv(mapCenter) {
var p3 = new google.maps.LatLng(mapCenter.lat(), mapCenter.lng()-0.5);
var p4 = new google.maps.LatLng(mapCenter.lat(), mapCenter.lng()+0.5);
return dist(p3,p4);
}
function dist(pt1, pt2) {
var dLat = ( pt2.lat() - pt1.lat() ) * Math.PI / 180;
var dLng = ( pt2.lng() - pt1.lng() ) * Math.PI / 180;
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(rad(pt1.lat())) * Math.cos(rad(pt2.lat())) *
Math.sin(dLng/2) * Math.sin(dLng/2);
var R = 6372800; //earth's radius
var distance = R * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
return distance;
}
function Area(polygon) {
var xPts=[];
for (i=1; i<polygon.length; i++) {
xPts[i-1] = ( polygon[i].lat() - polygon[0].lat() ) * latDiv;
}
var yPts=[];
for (i=1; i<polygon.length; i++) {
yPts[i-1] = ( polygon[i].lng() - polygon[0].lng() ) * lngDiv;
}
var area = 0;
j = polygon.length-2;
for (i=0; i<polygon.length-1; i++) {
area = area + ( xPts[j] + xPts[i] ) * ( yPts[j] - yPts[i] );
j = i;
}
area = Math.abs(area/2);
return area;
}
评论
您在此处尝试应用的多边形区域的笛卡尔公式不适用于在曲面(例如球体)上计算的距离和角度。通过使用纬度和经度,就好像它们是笛卡尔坐标一样,此公式也会产生另一个错误。可以考虑使用它的唯一情况就是(对于非常小的多边形而言)恰好是不必要使用hasrsine公式的情况。总的来说,此代码工作太辛苦了,没有收益。
– hu
2014年6月28日14:27
评论
您是否正在寻找一种可以做到这一点的工具(例如Esri的pe.dll)或公式?抱歉,我不确定...我在寻找公式。我将更新问题以更具体。
大量的数学示例位于此处计算纬度/经度点之间的距离,方位角以及更多信息,其中包括“目标位置”的解决方案给定点距起点的距离和方位”。
密切相关:gis.stackexchange.com/questions/2951/…。
这是链接到经纬度计算的页面[经纬度计算](movable-type.co.uk/scripts/latlong.html)也是此页面经纬度计算有一个代码+计算器