Arcpy根据点的访问和距已知点的距离计算坐标值

标签： Arcpy 分类： Gis 创建时间：2019-10-18 02:37:44 更新时间：2025-01-20 09:45:23

现在有一个需求，就是通过已知点的经纬度坐标，和另一个点相对于已知点的方向和距离，计算另一个点的经纬度坐标。好像网上的大部分都是根据两个经纬度坐标。

1.根据经纬度计算两点距离

根据Haversine公式，计算两者之间的距离。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#!/usr/bin/python
# -*- coding: UTF-8 -*-

import math
 
## km为单位
def distance(origin, destination):
    lat1, lon1 = origin
    lat2, lon2 = destination
    radius = 6371 # km
 
    dlat = math.radians(lat2-lat1)
    dlon = math.radians(lon2-lon1)
    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = radius * c
 
    return d

from math import sin, asin, cos, radians, fabs, sqrt

## m为单位
def haversine(lon1, lat1, lon2, lat2): # 经度1，纬度1，经度2，纬度2 （十进制度数）
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # 将十进制度数转化为弧度
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
 
    # haversine公式
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # 地球平均半径，单位为公里
    return c * r * 1000



if __name__=='__main__':
    # 38.35252313201552 (km为单位)
    print distance((22.599578, 113.973129), (22.6986848, 114.3311032))

    print distance((22.599578, 113.973129), (22.6986848, 114.3311032))
参考文章：
1.arcgis根据两个点经纬度，计算两点之间的距离
2.根据经纬度坐标计算距离-python
3.计算两个经纬度之间的距离(python算法)
4.PYTHON中如何计算2个坐标之间的距离

2.根据已知方向和距离计算另一个点的经纬度坐标

根据参考资料，我将c#代码转换成了python代码，虽然不懂原理，但是能运行啊。具体运行的效果，和真实值是有一定的误差的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#!/usr/bin/python
# -*- coding: UTF-8 -*-

import math

## 大地坐标系资料WGS-84 长半径a=6378137 短半径b=6356752.3142 扁率f=1/298.2572236

## 长半径a=6378137 
a = 6378137
## 短半径b=6356752.3142
b = 6356752.3142
## 扁率f=1/298.2572236
f = 1 / 298.2572236

## lon 经度，lat 纬度，brng 方位角，dist 距离(m)
def computerThatLonLat(lon,lat,brng,dist):
    alpha1 = rad(brng)
    sinAlpha1 = math.sin(alpha1)
    cosAlpha1 = math.cos(alpha1)

    tanU1 = (1 - f) * math.tan(rad(lat))
    cosU1 = 1 / math.sqrt((1 + tanU1 * tanU1))
    sinU1 = tanU1 * cosU1
    sigma1 = math.atan2(tanU1, cosAlpha1)
    sinAlpha = cosU1 * sinAlpha1
    cosSqAlpha = 1 - sinAlpha * sinAlpha
    uSq = cosSqAlpha * (a * a - b * b) / (b * b)
    A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))

    cos2SigmaM=0
    sinSigma=0
    cosSigma=0
    sigma = dist / (b * A)
    sigmaP = 2 * math.pi
    while math.fabs(sigma - sigmaP) > 1e-12:
        cos2SigmaM = math.cos(2 * sigma1 + sigma)
        sinSigma = math.sin(sigma)
        cosSigma = math.cos(sigma)
        deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)
                - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)))
        sigmaP = sigma
        sigma = dist / (b * A) + deltaSigma
    
    tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1
    lat2 = math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1, (1 - f) * math.sqrt(sinAlpha * sinAlpha + tmp * tmp))
    lamb = math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1)
    C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
    L = lamb - (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))

    lon=lon+deg(L)
    lat=deg(lat2)

    print lon
    print lat

## 度换成弧度,d 弧度
def rad(d):
    return d * math.pi / 180.0
## 弧度换成度,x 弧度
def deg(x):
    return x * 180 / math.pi

if __name__=='__main__':

    lon=120.353842
    lat=30.338461
    brng=90
    dist=1

    computerThatLonLat(lon, lat, brng, dist)
参考文章：
1.经纬度 距离，角度 计算另一点经纬度
2.已知经纬度，方位角，距离，求另一点经纬度
**本文已结束，感谢您的观看**
鼓励一下

往期推荐

文章目录
  1. 1. 1.根据经纬度计算两点距离
  2. 2. 2.根据已知方向和距离计算另一个点的经纬度坐标
小额赞助
本人提供免费与付费咨询服务，感谢您的支持！赞助请发邮件通知，方便公布您的善意！
**光 3.01 元
Sun 3.00 元
bibichuan 3.00 元
捐赠者名单 合计：9.01 元
微信公众号
广告位
诚心邀请广大金主爸爸洽谈合作
每日一省
isNaN 和 Number.isNaN 函数的区别？

1.函数 isNaN 接收参数后，会尝试将这个参数转换为数值，任何不能被转换为数值的的值都会返回 true，因此非数字值传入也会返回 true ，会影响 NaN 的判断。

2.函数 Number.isNaN 会首先判断传入参数是否为数字，如果是数字再继续判断是否为 NaN ，不会进行数据类型的转换，这种方法对于 NaN 的判断更为准确。

2021-04-07 链接地址
每日二省
为什么0.1+0.2 ! == 0.3，如何让其相等?

一个直接的解决方法就是设置一个误差范围，通常称为“机器精度”。对JavaScript来说，这个值通常为2-52，在ES6中，提供了Number.EPSILON属性，而它的值就是2-52，只要判断0.1+0.2-0.3是否小于Number.EPSILON，如果小于，就可以判断为0.1+0.2 ===0.3。

2021-04-07 链接地址
每日三省
== 操作符的强制类型转换规则？

1.首先会判断两者类型是否**相同，**相同的话就比较两者的大小。

2.类型不相同的话，就会进行类型转换。

3.会先判断是否在对比 null 和 undefined，是的话就会返回 true。

4.判断两者类型是否为 string 和 number，是的话就会将字符串转换为 number。

5.判断其中一方是否为 boolean，是的话就会把 boolean 转为 number 再进行判断。

6.判断其中一方是否为 object 且另一方为 string、number 或者 symbol，是的话就会把 object 转为原始类型再进行判断。

2021-04-07 链接地址
每日英语
Happiness is time precipitation, smile is the lonely sad.
幸福是年华的沉淀，微笑是寂寞的悲伤。
2016 ——— 美文