首页 > 解决方案 > How to check if a location is inside a circle on a map

问题描述

I'm trying to check if a location lat/long is inside a circle on a map using Matlab. I used this formula :

d=R*acos(sin(S_lat(t)*d2r)*sin(N_lat*d2r)+cos(S_lat(t)*d2r)*cos(N_lat*d2r)*cos((S_long(t)-N_long)*d2r))

and the following formula:

a = sin((S_lat(t)-N_lat)*d2r / 2)^2 + cos(N_lat*d2r) * cos(S_lat(t)*d2r) * sin((S_long(t)-N_long)*d2r / 2)^2;
c = 2 * asin(sqrt(a));
d = R * c;

But it didn't give me the right answers.

标签: matlabdictionarygeometry

解决方案


与此页面相比,您的第二个距离公式看起来几乎正确

a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
c = 2 ⋅ atan2( √a, √(1−a) )
d = R ⋅ c

据我所知,对于 0..1 范围内的正参数,应该给出与第二行中的公式asin相同的结果。atan2如果 的值a可能稍微超出这个区间,反正弦会给出错误的结果,所以

如果 atan2 不可用,则 c 可以从 2 ⋅ asin( min(1, √a) ) 计算(包括防止舍入错误)。

你能举一个有错误的数据点的例子吗?

Python代码(ideone玩)

import math
def eadist(lat1, lon1, lat2, lon2):
    print(lat1, lon1, lat2, lon2)
    d2r = math.pi / 180
    R = 6378
    dlat = (lat2 -lat1) * d2r
    dlon = (lon2 -lon1) * d2r
    a = (math.sin(dlat/2))**2 + math.cos(lat1*d2r)*math.cos(lat2*d2r) * (math.sin(dlon/2))**2
    print("a=",a)
    c = 2*math.asin(math.sqrt(a))
    print("c=",c)
    dist = c * R
    return dist

print(eadist(45, 90, 45, 91))
print(eadist(45, 90, 46, 90))
print(eadist(45, 90, 46, 91))
print(eadist(45, 90, -45, -90))

给出正确的结果

45 90 45 91
a= 3.8076210902190215e-05
c= 0.012341263173265506
78.71257651908739  //1 degree by parallel

45 90 46 90
a= 7.615242180438042e-05
c= 0.017453292519943295
111.31709969219834   //1 degree by meridian

45 90 46 91
a= 0.0001135583120069672
c= 0.021313151879913415
135.93528269008777   //diagonal step

45 90 -45 -90
a= 1.0
c= 3.141592653589793
20037.0779445957   //antipodal point, a half of meridian length

推荐阅读