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

第二种方法：

 1 #import <Foundation/Foundation.h>
 2 
 3 @interface MyLatLng : NSObject {
 4     double m_LoDeg,m_LoMin,m_LoSec;
 5     double m_LaDeg,m_LaMin,m_LaSec;
 6     double m_Longitude,m_Latitude;
 7     double m_RadLo,m_RadLa;
 8     double Ec;
 9     double Ed;
10 }
11 - (id)init:(double)longitude latitude:(double)latitude;
12 
13 @property (assign, nonatomic) double m_LoDeg;
14 @property (assign, nonatomic) double m_LoMin;
15 @property (assign, nonatomic) double m_LoSec;
16 @property (assign, nonatomic) double m_LaDeg;
17 @property (assign, nonatomic) double m_LaMin;
18 @property (assign, nonatomic) double m_LaSec;
19 @property (assign, nonatomic) double m_Longitude;
20 @property (assign, nonatomic) double m_Latitude;
21 @property (assign, nonatomic) double m_RadLo;
22 @property (assign, nonatomic) double m_RadLa;
23 @property (assign, nonatomic) double Ec;
24 @property (assign, nonatomic) double Ed;
25 @end

 1 #import "MyLatLng.h"
 2 #define RC 6378137
 3 #define RJ 6356725
 4 
 5 @implementation MyLatLng
 6 @synthesize m_LoDeg;
 7 @synthesize m_LoMin;
 8 @synthesize m_LoSec;
 9 @synthesize m_LaDeg;
10 @synthesize m_LaMin;
11 @synthesize m_LaSec;
12 @synthesize m_Longitude;
13 @synthesize m_Latitude;
14 @synthesize m_RadLa;
15 @synthesize m_RadLo;
16 @synthesize Ec;
17 @synthesize Ed;
18 
19 - (id)init:(double)longitude latitude:(double)latitude{
20     self = [super init];
21     
22     if (self) {
23         m_LoDeg=(int)longitude;
24         m_LoMin=(int)((longitude-m_LoDeg)*60);
25         m_LoSec=(longitude-m_LoDeg-m_LoMin/60.)*3600;
26         
27         m_LaDeg=(int)latitude;
28         m_LaMin=(int)((latitude-m_LaDeg)*60);
29         m_LaSec=(latitude-m_LaDeg-m_LaMin/60.)*3600;
30         
31         m_Longitude=longitude;
32         m_Latitude=latitude;
33         m_RadLo=longitude*M_PI/180.;
34         m_RadLa=latitude*M_PI/180.;
35         Ec=RJ+(RC-RJ)*(90.-m_Latitude)/90.;
36         Ed=Ec*cos(m_RadLa);
37     }
38     
39     return self;
40 }
41 @end

 1 #import "ViewController.h"
 2 #import "MyLatLng.h"
 3 
 4 @interface ViewController ()
 5 
 6 @end
 7 
 8 @implementation ViewController
 9 
10 - (void)viewDidLoad {
11     [super viewDidLoad];
12     // Do any additional setup after loading the view, typically from a nib.
13     
14     MyLatLng *A = [[MyLatLng alloc]init:113.249648 latitude:23.401553];
15     MyLatLng *B = [[MyLatLng alloc]init:113.246033 latitude:23.403362];
16     
17     [self getAngle:A B:B];
18 }
19 
20 - (double)getAngle:(MyLatLng *)A B:(MyLatLng *)B {
21     double dx=(B.m_RadLo-A.m_RadLo)*A.Ed;
22     double dy=(B.m_RadLa-A.m_RadLa)*A.Ec;
23     double angle=0.0;
24     angle = atan(fabs(dx/dy))*180./M_PI;
25     double dLo=B.m_Longitude-A.m_Longitude;
26     double dLa=B.m_Latitude-A.m_Latitude;
27     if(dLo>0&&dLa<=0){
28         angle=(90.-angle)+90;
29     }
30     else if(dLo<=0&&dLa<0){
31         angle=angle+180.;
32     }else if(dLo<0&&dLa>=0){
33         angle= (90.-angle)+270;
34     }
35     return angle;
36 }
37 
38 - (void)didReceiveMemoryWarning {
39     [super didReceiveMemoryWarning];
40     // Dispose of any resources that can be recreated.
41 }