我正在使用Matlab检查一个位置lat/long是否在地图上的圆圈内。我用了这个公式：
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;
但它并没有给我正确的答案。
R=6378；%地球引力d2r=π/180；%角度到弧度
N_ lat=45；N_ long=90；
rov=70；
[latc，longc]=scircle1(N_lat，N_long，rov(；%取一个以(N_lat，N_long(为中心的圆，并按其半径
D=零(1，长度(latc((；%检查从圆心到其任意点的距离%周长等于rov使用这个公式
对于i=1：长度(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));
%公式2
%D(i(=R*sqrt((latc(i(-N-lat(^2+(longc(i，-N_long(^2(；
%公式3%a=sin((latc(i(-N-lat(d2r/2(^2+cos(N_latd2r(*cos(latc；%c=2*asin(sqrt(a((；%D(i(=R*c；%公式4%D(i(=R距离(latc(i(，longc(i，N_lat，N_long(；
结束
按矢量几何和球坐标：
您可以通过将经度/纬度角度转换为单位球体上的笛卡尔坐标
X = cos(Θ) sin(φ), Y = sin(Θ) sin(φ), Z = cos(φ)
您可以为测试点和圆心计算这些值。然后点积的弧余弦给出了中心角，你可以将其与地球半径上的圆半径进行比较。
*圆弧长度，而非圆平面内的半径。
与本页相比，您的第二个距离公式看起来几乎正确
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代码(可播放的视频(
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
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,N_long,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