python之由方向和距离求坐标(二分逼近法)

15195775117

问题由来：从已知坐标往X方向L米处的坐标是多少？
目前，我只知道距离怎么求，求坐标没有现成函数，于是我利用二分法来逼近精确坐标。

思路是，我先选坐标原点和大约1万米外的点(近似为经纬度0.1)为直线两端，
然后选择端点的中间点为新端点，
然后新端点和旧端点继续折中出新端点，......，
直到其中一个端点与原点的距离在L±1m的精度范围内。


from numpy import sin,cos
from geopy.distance import geodesic
from numpy import pi


#该函数用于求从某点开始，朝某方向一定距离的位置
#例如从[120.12813,30.28718]往东北方向500米的位置：
# theta=45，东北方向，笛卡尔坐标系角度，角度制，从正东逆时针旋转
# core=[120.12813,30.28718]#坐标
# distance=500 #500米远


def neighbor(core,theta,distance):
    theta2=theta/180.0*pi  #转换为弧度制
    vec=[cos(theta2),sin(theta2)]  #vec是单位矢量，代表方向
   
    constraint=1#精确度，单位米，不达到这个精度循环不停止
   
    x1=float(core[0])  #原点经度
    y1=float(core[1])  #原点纬度
    x2=x1+0.1*vec[0]  #从约1万米远为初始“远端点”
    y2=y1+0.1*vec[1]
    mx=(x1+x2)/2  #初始“折中点”
    my=(y1+y2)/2
   
    #远近2个端点与原点的距离：
    d1=geodesic((y1,x1),(core[1],core[0])).m  #注意，该函数先纬度后经度！
    d2=geodesic((y2,x2),(core[1],core[0])).m
   
    #循环，直到端点之一的距离符合精度：
    while abs(d1-distance) > constraint and abs(d2-distance) > constraint:
        #如果折中点的距离过远，则折中点为新的远端点：
        if geodesic((my,mx),(core[1],core[0])).m > distance+constraint:
            x2=mx
            y2=my
            mx=(x1+x2)/2
            my=(y1+y2)/2
   
        #如果折中点的距离过远，则折中点为新的近端点：
        if geodesic((my,mx),(core[1],core[0])).m < distance+constraint:
            x1=mx
            y1=my
            mx=(x1+x2)/2
            my=(y1+y2)/2
        
        #结束一轮端点重设后，计算2个端点与原点的距离，判断循环是否结束：
        d1=geodesic((y1,x1),(core[1],core[0])).m
        d2=geodesic((y2,x2),(core[1],core[0])).m
   
    #三元表达式，选择(x1,y1)和(x2,y2)两点中最接近距离要求的
    x=x1 if d1 < d2 else x2
    y=y1 if d1 < d2 else y2
   
    err=geodesic((y,x),(core[1],core[0])).m-distance  #距离误差
    if err > constraint:
        print('误差大于约束值，出错！')
    # print('坐标=',x,y)
    # print('误差=','%.2f'%err,'m')
    return [x,y,err]


if __name__ == "__main__":
    #现在绘制22.5度×100m的网格
    import numpy as np
    import matplotlib.pyplot as plt
    theta=np.arange(0,361,22.5)  #角度序列
    distance=np.arange(100,1000,100)  #距离序列
    n=len(theta)*len(distance)  #求个数
    theta,distance=np.meshgrid(theta,distance)  #网格化
    theta=theta.reshape(n,1)  #转为1维
    distance=distance.reshape(n,1)  #转为1维
    theta_dis=np.hstack((theta,distance))  #并排堆积，用于后面的循环
    core=[120.12813,30.28718]  #圆心
    book=[0,0,0]  #存储结果的变量
    for [i,j] in theta_dis:
        result=neighbor(core,i,j)
        book=np.vstack((book,result))
    book=book[1:,:]  #舍去初始值
    plt.axes([0,0,1,1])  #撑满画布
    plt.scatter(core[0],core[1],marker='*')  #绘制圆心
    plt.scatter(book[:,0],book[:,1])  #绘制散点
    plt.axis('off')  #隐藏坐标轴
    plt.show()

15195775117

该算法的IDL版见《IDL之由方向和距离求坐标(二分逼近法)》：
http://bbs.06climate.com/forum.p ... p;extra=#pid1079766

15195775117

周末浏览了高德地图开发文档，它也提供了类似的功能：
根据距离差计算另一个经纬度：
var lnglat3 = lnglat1.offset(100,50)//lnglat1向东100m，向北50m的位置的经纬度
来源在此