iOS 根据经纬度计算与地理北极夹角

http://www.aiuxian.com/article/p-2767848.html

#define toDeg(X) (X*180.0/M_PI)

/**

 * @method 根据两点经纬度，计算与真北方夹角

 *

 * @param longitude1

 * @param latitude1

 * @param longitude2

 * @param latitude2 // 目标点

 */

- (void)getAngle:(double)longitude1 latitude1:(double)latitude1 longitude2:(double)longitude2 latitude2:(double)latitude2 {

    double cos_c = cos(90 - latitude2)*cos(90 - latitude1) + sin(90 - latitude2)*sin(90 - latitude1)*cos(longitude2-longitude1);

    double sin_c = sqrt(1 - pow(cos_c, 2));

    double z = asin(sin(90 - latitude2)*sin(longitude2 - longitude1)/sin_c);

    z = toDeg(z);

// A（起始点）为原点B目标点

    if(longitude1 < longitude2 && latitude1 < latitude2) { // 第一象限

    } else if (longitude1 < longitude2 && latitude1 > latitude2){ // 第二象限

        z += 360;

    } else { // 三四象限

        z = 180 - z;

    }

    cityHeading = z;

    NSLog(@"城市夹角:%f",z);

}

此处设定求B相对于A的方位角，即A为当前位置，B为目标位置

Aj：A点经度

Aw：A点纬度

Bj：B点经度

Bw：B点纬度

北纬为正，南纬为负；东经为正，西经为负

经纬度使用度，DDD.DDDDDD°，非度分或度分秒。

度数未加说明均采用角度制

R：地球平均半径

Azimuth：方位角，以真北为0度起点，由东向南向西顺时针旋转360度

这里需要注意一点，我们一开始的假设便是求B点相对于A点的方位角，因此这里是Bj-Aj，不要写反，否则得不到正确结果。

算到这里，还没有完，得到的结果并不总符合我们对方位角的定义，因此要根据B相对于A的位置在四个象限两个轴上进行讨论，依据不同情况对计算结果进行不同处理。假设A点固定于原点，则:

B点在第一象限，Azimuth=A；

B在第二象限，Azimuth=360+A;

B在第三四象限，Azimuth=180-A。

这里只说了象限的讨论结果，因为轴上的讨论更复杂些，要结合程序运行环境一起考虑，考虑的主要因素是系统的计算精度。譬如，在三面角余弦公式中，当AB点纬度值相同时，对公式的值起决定作用的就是cos(Bj-Aj)这一项，当Bj-Aj的值比较小时，例如0.0001（这在赤道地区对应的长度为11米左右），用一般的计算器计算时值为1，这样，后面的计算便不可能完成。但是，如果用计算机计算则为0.999999999998476913…………。所以，基于以上原因，需要对轴的“范围进项扩充”，要用单片机、手机运算的尤其要注意。

经过一系列计算，最后，就得到了最终结果。

似乎有人注意到了，以上的计算都是把地球看成标准的球体，而事实是地球是个椭圆，其实，地球的偏心率极低，各位可以将此法得到的计算结果与谷歌地球（WGS84坐标系统，我说的不是谷歌地图）上的结果进行对比，偏差是非常小的（我测的几个值，最大偏差0.5度）。

二、距离的求算

其实，“眼尖”的或许已经注意到了，第一步的余弦值结果就可以直接用来求算AB两点间的球面距离，用反余弦函数求得c的度数，再将度数转换为弧度，乘以地球半径就得到了两点间的球面距离。

公式为

这里要注意，L的单位与R的单位一致，单位不同的不要忘记换算

短距离（例如100米，30米）使用这个公式，计算出的结果与谷歌地球给出的距离偏差在0.5%以下，长距离计算时，偏差则可以降至0.01%以下。求算的距离越大，偏差越小，就是这个公式的特点，原因不说自明。

PS:对于一些GPS接收机，其数据格式为NMEA-0183，经纬度数据为DDDMM.MMMM，需要将它转换为度，公式为：

经纬度（度）=DDD+MM.MMMM/60

第二种方法：

#import <Foundation/Foundation.h>

@interface MyLatLng : NSObject {
double m_LoDeg,m_LoMin,m_LoSec;
double m_LaDeg,m_LaMin,m_LaSec;
double m_Longitude,m_Latitude;
double m_RadLo,m_RadLa;
double Ec;
double Ed;
}
- (id)init:(double)longitude latitude:(double)latitude;

@property (assign, nonatomic) double m_LoDeg;
@property (assign, nonatomic) double m_LoMin;
@property (assign, nonatomic) double m_LoSec;
@property (assign, nonatomic) double m_LaDeg;
@property (assign, nonatomic) double m_LaMin;
@property (assign, nonatomic) double m_LaSec;
@property (assign, nonatomic) double m_Longitude;
@property (assign, nonatomic) double m_Latitude;
@property (assign, nonatomic) double m_RadLo;
@property (assign, nonatomic) double m_RadLa;
@property (assign, nonatomic) double Ec;
@property (assign, nonatomic) double Ed;
@end

#import "MyLatLng.h"
#define RC 6378137
#define RJ 6356725

@implementation MyLatLng
@synthesize m_LoDeg;
@synthesize m_LoMin;
@synthesize m_LoSec;
@synthesize m_LaDeg;
@synthesize m_LaMin;
@synthesize m_LaSec;
@synthesize m_Longitude;
@synthesize m_Latitude;
@synthesize m_RadLa;
@synthesize m_RadLo;
@synthesize Ec;
@synthesize Ed;

- (id)init:(double)longitude latitude:(double)latitude{
self = [super init];

if (self) {
m_LoDeg=(int)longitude;
m_LoMin=(int)((longitude-m_LoDeg)*60);
m_LoSec=(longitude-m_LoDeg-m_LoMin/60.)*3600;

m_LaDeg=(int)latitude;
m_LaMin=(int)((latitude-m_LaDeg)*60);
m_LaSec=(latitude-m_LaDeg-m_LaMin/60.)*3600;

m_Longitude=longitude;
m_Latitude=latitude;
m_RadLo=longitude*M_PI/180.;
m_RadLa=latitude*M_PI/180.;
Ec=RJ+(RC-RJ)*(90.-m_Latitude)/90.;
Ed=Ec*cos(m_RadLa);
}

return self;
}
@end

#import "ViewController.h"
#import "MyLatLng.h"

@interface ViewController ()

@end

@implementation ViewController

- (void)viewDidLoad {
[super viewDidLoad];
// Do any additional setup after loading the view, typically from a nib.

MyLatLng *A = [[MyLatLng alloc]init:113.249648 latitude:23.401553];
MyLatLng *B = [[MyLatLng alloc]init:113.246033 latitude:23.403362];

[self getAngle:A B:B];
}

- (double)getAngle:(MyLatLng *)A B:(MyLatLng *)B {
double dx=(B.m_RadLo-A.m_RadLo)*A.Ed;
double dy=(B.m_RadLa-A.m_RadLa)*A.Ec;
double angle=0.0;
angle = atan(fabs(dx/dy))*180./M_PI;
double dLo=B.m_Longitude-A.m_Longitude;
double dLa=B.m_Latitude-A.m_Latitude;
if(dLo>0&&dLa<=0){
angle=(90.-angle)+90;
}
else if(dLo<=0&&dLa<0){
angle=angle+180.;
}else if(dLo<0&&dLa>=0){
angle= (90.-angle)+270;
}
return angle;
}

- (void)didReceiveMemoryWarning {
[super didReceiveMemoryWarning];
// Dispose of any resources that can be recreated.
}