根据经纬度计算两点之间的距离

#define PI 3.1415926
double LantitudeLongitudeDist(double lon1,double lat1,
double lon2,double lat2)
{
double er = 6378137; // 6378700.0f;
//ave. radius = 6371.315 (someone said more accurate is 6366.707)
//equatorial radius = 6378.388
//nautical mile = 1.15078
double radlat1 = PI*lat1/180.0f;
double radlat2 = PI*lat2/180.0f;
//now long.
double radlong1 = PI*lon1/180.0f;
double radlong2 = PI*lon2/180.0f;
if( radlat1 < 0 ) radlat1 = PI/2 + fabs(radlat1);// south
if( radlat1 > 0 ) radlat1 = PI/2 - fabs(radlat1);// north
if( radlong1 < 0 ) radlong1 = PI*2 - fabs(radlong1);//west
if( radlat2 < 0 ) radlat2 = PI/2 + fabs(radlat2);// south
if( radlat2 > 0 ) radlat2 = PI/2 - fabs(radlat2);// north
if( radlong2 < 0 ) radlong2 = PI*2 - fabs(radlong2);// west
//spherical coordinates x=r*cos(ag)sin(at), y=r*sin(ag)*sin(at), z=r*cos(at)
//zero ag is up so reverse lat
double x1 = er * cos(radlong1) * sin(radlat1);
double y1 = er * sin(radlong1) * sin(radlat1);
double z1 = er * cos(radlat1);
double x2 = er * cos(radlong2) * sin(radlat2);
double y2 = er * sin(radlong2) * sin(radlat2);
double z2 = er * cos(radlat2);
double d = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
//side, side, side, law of cosines and arccos
double theta = acos((er*er+er*er-d*d)/(2*er*er));
double dist = theta*er;
return dist;
}
### 基于经纬度计算两点距离的编程方法 #### 背景介绍 在地理信息系统(GIS)以及空间数据分析领域,经常需要根据地球表面两点经纬度计算它们之间的球面距离。这种需求可以通过多种方式实现,其中最常用的是 **Haversine Formula** 和 **Vincenty's Formulae**。 #### Haversine 公式的理论基础 Haversine 公式是一种用于计算地球上任意两点之间大圆距离的有效算法[^1]。该公式的具体表达形式如下: \[ a = \sin^2\left(\frac{\Delta\phi}{2}\right) + \cos(\phi_1)\cdot\cos(\phi_2)\cdot\sin^2\left(\frac{\Delta\lambda}{2}\right) \] \[ c = 2\cdot\textrm{atan2}(\sqrt{a}, \sqrt{1-a}) \] \[ d = R\cdot c \] 其中: - \(R\) 是地球半径(通常取平均值约为 6371 km) - \(\phi_1, \phi_2\) 表示两点的纬度(单位为弧度) - \(\Delta\phi\) 表示纬度差 - \(\Delta\lambda\) 表示经度差 需要注意的是,在实际应用中,角度需转换为弧度制才能代入上述公式进行运算。 #### MATLAB 实现代码 以下是基于 Haversine 公式的 MATLAB 实现代码: ```matlab function distance = haversine(lat1, lon1, lat2, lon2) % 将角度转化为弧度 radLat1 = deg2rad(lat1); radLon1 = deg2rad(lon1); radLat2 = deg2rad(lat2); radLon2 = deg2rad(lon2); % 地球半径 (km) R = 6371; % 计算纬度和经度差异 deltaLat = radLat2 - radLat1; deltaLon = radLon2 - radLon1; % 应用 Haversine 公式 a = sin(deltaLat / 2)^2 + cos(radLat1) * cos(radLat2) * sin(deltaLon / 2)^2; c = 2 * atan2(sqrt(a), sqrt(1 - a)); distance = R * c; % 返回距离 (km) end ``` 此函数接受四个参数:起点和终点的纬度与经度,并返回两者之间的球面距离(单位为公里)。注意输入的角度应为十进制度数而非其他格式。 #### Python 实现代码 Python 中可以借助 `geopy` 或者手动编写 Haversine 函数完成相同功能。下面展示一种不依赖外部库的手动实现方案: ```python import math def calculate_distance_haversine(lat1, lon1, lat2, lon2): # 地球半径 (km) R = 6371 # 将角度转为弧度 phi1 = math.radians(lat1) phi2 = math.radians(lat2) delta_phi = math.radians(lat2 - lat1) delta_lambda = math.radians(lon2 - lon1) # 使用 Haversine 公式 a = math.sin(delta_phi / 2)**2 + \ math.cos(phi1) * math.cos(phi2) * math.sin(delta_lambda / 2)**2 c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a)) distance = R * c # 结果以千米表示 return round(distance, 2) # 保留两位小数并返回 ``` 如果希望简化开发过程,则可采用第三方库如 geopy 提供的功能[^2]: ```python from geopy.distance import geodesic point1 = (lat1, lon1) # 替换为实际坐标 point2 = (lat2, lon2) # 替换为实际坐标 distance_km = geodesic(point1, point2).kilometers print(f"Distance: {round(distance_km, 2)} km") ``` 这种方法更加简洁明了,适合快速原型设计阶段使用。 ---
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值