objc

[Objective-C] 特定の緯度経度と角度と距離を元に緯度経度を算出する(楕円あり)

サイトウ

2015年11月14日 • 2 min read

現在地を起点として、特定距離の矩形を作成する。
みたいなことが必要だったので計算式を書いた(ほぼコピペ)。。。けど使われないことになった。

やってることとしては

ある地点を中心とした半径Xメートル同心円上における角度Aの地点の緯度経度を求めよ(ただし地球は楕円形として扱うこと)

みたいな感じでしょうか。文系な自分にはこういうのよくわからないです。。

結果として下記ソースでほぼイケるけど、浮動小数点数の演算ばかりなので誤差が気になるところ。メソッドのコメントに記載しているURLでJavascriptのソースがあったので、Objective-Cに書き直しただけって感じですが、一部変えています。詳細はソース見てもらえれば。

参考サイトではdo whileでループしているが

この部分の動作を見てみたら、Javascriptでは精度を上げるようなループをしていたけど期待通りに実装する方法がわからなかったので、ループの中身を1度実行するだけに変更しています。

#import <CoreLocation/CoreLocation.h>

#define EARTH           6378137.f
#define EARTH_FLAT      1.f / 298.257222101
#define EARTH_SHORT     EARTH * (1.f - EARTH_FLAT)

#define RADIANS_TO_DEGREES(radians) (radians * 180.f / M_PI)
#define DEGREES_TO_RADIANS(degrees) (degrees / 180.f * M_PI)

@implementation AppDelegate

/**
 @param point 起点となる緯度経度
 @param angle 角度、北を0度として時計回りに360度
 @param distance 距離(メートル)
 */
+ (CLLocationCoordinate2D)calcLocationWithPoint:(CLLocationCoordinate2D)point
                                          angle:(double)angle
                                       distance:(double)distance {
    /*
     http://tancro.e-central.tv/grandmaster/script/vincentyJS.html
     http://tancro.e-central.tv/grandmaster/excel/vincenty.html
     
     MKCoordinateRegionが使えるかと思ったけど
     MKMapRectが無いと期待している値が取れなかった
     */
    // 現在地の経度
    double latitude = DEGREES_TO_RADIANS(point.latitude);
    // 現在地の緯度
    double longitude = DEGREES_TO_RADIANS(point.longitude);
    // 算出する地点の角度
    double alpha12 = DEGREES_TO_RADIANS(angle);
    
    double U1 = atan((1 - EARTH_FLAT) * tan(latitude));
    double sigma1 = atan(tan(U1) / cos(alpha12));
    double alpha = asin(cos(U1) * sin(alpha12));
    double u2 = pow(cos(alpha), 2) * (pow(EARTH, 2) - pow(EARTH_SHORT, 2)) / pow(EARTH_SHORT, 2);
    double A = 1 + (u2 / 16384.f) * (4096.f + u2 * (-768 + u2 * (320.f - 175.f * u2)));
    double B = (u2 / 1024.f) * (256.f + u2 * (-128 + u2 * (74.f - 47.f * u2)));
    double sigma = distance / EARTH_SHORT / A;

    /*
     参考サイトではdo whileでループしているが
     iPhoneで試したところうまく動かなかったためループさせないよう修正
     */
    double dm2;
    dm2 = 2.f * sigma1 + sigma;
    double x1 = cos(sigma) * ( -1 + 2.f * pow(cos(dm2),2) ) - B / 6.f * cos(dm2) * ( -3.f + 4.f * pow(sin(dm2),2)) * ( -3.f + 4.f * pow(cos(dm2), 2));
    double dSigma = B * sin(sigma) * ( cos(dm2) + B / 4.f * x1);
    sigma = distance / EARTH_SHORT / A + dSigma;
    
    double x = sin(U1) * cos(sigma) + cos(U1) * sin(sigma) * cos(alpha12);
    double y = (1.f - EARTH_FLAT) * pow( pow( sin(alpha), 2) + pow(sin(U1) * sin(sigma) - cos(U1) * cos(sigma) * cos(alpha12), 2), .5);
    double lambda = atan(sin(sigma) * sin(alpha12) / (cos(U1) * cos(sigma) - sin(U1) * sin(sigma) * cos(alpha12)));
    double C = (EARTH_FLAT / 16.f) * pow(cos(alpha12), 2) * (4.f + EARTH_FLAT * (4.f - 3.f * pow(cos(alpha), 2)));
    double z = cos(dm2) + C * cos(sigma) * (-1 + 2.f * pow(cos(dm2), 2));
    double omega = lambda - (1.f - C) * EARTH_FLAT * sin(alpha) * (sigma + C * sin(sigma) * z);
    double lat = atan(x / y);
    double lng = longitude + omega;
    CLLocationCoordinate2D coordinate = CLLocationCoordinate2DMake(RADIANS_TO_DEGREES(lat),
                                                                   RADIANS_TO_DEGREES(lng));
    return coordinate;
}

# pragma mark - AppDelegate methods

- (BOOL)application:(UIApplication *)application didFinishLaunchingWithOptions:(NSDictionary *)launchOptions {
    // 東京駅付近の緯度経度
    CLLocationCoordinate2D basePoint = CLLocationCoordinate2DMake(35.681330, 139.766029);
    CLLocationCoordinate2D result = [[self class] calcLocationWithPoint:basePoint
                                                          angle:0.f
                                                       distance:100.f];
    NSLog(@"basePoint:  %f, %f", basePoint.latitude, basePoint.longitude);// basePoint:  35.681330, 139.766029
    NSLog(@"result:     %f, %f", result.latitude, result.longitude); // result:     35.682225, 139.766029

    return YES;
}


@end

結果

iOSでこういった計算を行うとき、MapKit.frameworkのMKMapViewを使用しているとMKRectが使えるのでラクになるけど、MapViewを作ってるのはちょっとなぁと思ったので計算式を書きましたという経緯です。

参考サイト

Sign up for more like this.