如何检查位置是否在地图上的圆圈内

如何解决如何检查位置是否在地图上的圆圈内

我正在尝试使用Matlab检查纬度/经度是否在地图上的圆内。我用了这个公式：

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))

和以下公式：

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;

但是它没有给我正确的答案。

解决方法

与this page相比，您的第二个距离公式看起来几乎是正确的

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

据我所知，对于范围为0..1的正参数，asin的结果应与第二行中的atan2公式相同。 如果a的值可能略微超出此间隔，则反正弦将给出错误的结果，所以

如果atan2不可用，则可以从2⋅asin（min（1， √a））（包括防止舍入错误）。

您能举例说明有错误的数据点吗？

Python代码（ideone to play with)

import math
def eadist(lat1,lon1,lat2,lon2):
    print(lat1,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,46,90))
print(eadist(45,-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

,

R = 6378; ％地球 d2r = pi / 180; ％度到弧度

N_lat = 45; N_long = 90;

rov = 70;

[latc，longc] = scircle1（N_lat，N_long，rov）;％以（N_lat，N_long）为中心对一个圆进行半径补偿，并rov为其半径

D = zeros（1，length（latc））; ％检查从圆心到圆的任何点的距离 使用此公式，％perimeter等于rov

对于我= 1：length（latc）

% Formula1
%D(i)=acos(sin(latc(i)*d2r)*sin(N_lat*d2r)+cos(latc(i)*d2r)*cos(N_lat*d2r)*cos((longc(i)-N_long)*d2r));

％Formula2
％D（i）= R * sqrt（（latc（i）-N_lat）^ 2 +（longc（i）-N_long）^ 2）;

％Formula3 ％a = sin（（latc（i）-N_lat） d2r / 2）^ 2 + cos（N_lat d2r）* cos（latc（i）* d2r）* sin（（longc（i） -N_long） d2r / 2）^ 2; ％c = 2 * asin（sqrt（a））; ％D（i）= R * c; ％Formula4 ％D（i）= R 距离（latc（i），longc（i），N_lat，N_long）;

结束

,

通过矢量几何和球形坐标：

您可以通过以下方式将经度/纬度角转换为单位球面上的直角坐标：

X = cos(Θ) sin(φ),Y = sin(Θ) sin(φ),Z = cos(φ)

您可以为测试点和圆心计算这些值。然后，点积的反余弦值将为您提供圆心角，您可以将其与地球半径上的圆弧*进行比较。

---

*圆弧长度，而不是圆平面中的半径。

,

 R=6378;        %Earth raius
 d2r = pi/180;  %degrees to radians  

N_lat=45;
N_long=45;
rov=50;

%World map
        worldmap world
        load coastlines
        [latcells,loncells] = polysplit(coastlat,coastlon);
        plotm(coastlat,coastlon,'green')
   
hold on;
plotm(N_lat,N_long,'b*');
hold on;
[latc,longc] = scircle1(N_lat,rov);%compute a circle with (N_lat,N_long)as a center and rov its radius
plotm(latc,longc,'r-');

D=zeros(1,length(latc));
%Check if the distance from the center of the circle to any point of its
%perimeter is equal to rov using this formulas

for i =1:length(latc)
%     
%     % Formula1
%     %D(i)=acos(sin(latc(i)*d2r)*sin(N_lat*d2r)+cos(latc(i)*d2r)*cos(N_lat*d2r)*cos((longc(i)-N_long)*d2r));
%     
% %  Formula2        
% %  D(i)=R*sqrt((latc(i)-N_lat)^2+(longc(i)-N_long)^2);
% 
% %Formula3
% % a = sin((latc(i)-N_lat)*d2r / 2)^2 + cos(N_lat*d2r) * cos(latc(i)*d2r) * sin((longc(i)-N_long)*d2r / 2)^2;
% % c = 2 * asin(sqrt(a));
% % D(i) = R * c;
% 
% %Formula4
% 
% D(i)=distance('gc',latc(i),longc(i),N_lat,N_long);

% %Formula5
% x = (longc(i)-N_long)*d2r*cos((latc(i)+N_lat)*d2r/2);
% y = (latc(i)-N_lat)*d2r; 
% D(i)=R*sqrt(x^2+y^2);

%Formula6
 
x=cos(latc(i)*d2r)*cos(longc(i)*d2r)-cos(N_lat*d2r)*cos(N_long*d2r);
y=cos(latc(i)*d2r)*sin(longc(i)*d2r)-cos(N_lat*d2r)*sin(N_long*d2r);
z=sin(latc(i)*d2r)-sin(N_lat*d2r);

c=sqrt(x^2+y^2+z^2);
angle=asin(c/2);
D(i)=R*angle;

end

    enter code here