(long1, lat1), (long2, lat2)两点间的距离计算和中点计算

This article is categorized as "Garbage" . It should NEVER be appeared in your search engine's results.

Table of Contents

(long1, lat1), (long2, lat2)两点间的距离计算

之前一直用long, lat -> utm (x, y)的方法计算直线距离，但忽略了一个问题：utm可能跨区域。之前用到这个计算距离的场景都是小范围距离计算，但现在不能这么用了，现在肯定要用通用的标准方法。

标准算法是 Haversine Formula . 算法本身非常容易实现：

🔗 [engineering.purdue.edu/ece477/Archive/2008/Spring/S08-Grp03/nb/datasheets/haversine.pdf] https://engineering.purdue.edu/ece477/Archive/2008/Spring/S08-Grp03/nb/datasheets/haversine.pdf

https://engineering.purdue.edu/ece477/Archive/2008/Spring/S08-Grp03/nb/datasheets/haversine.pdf

所以自己就能写：

import math


# 注意返回的单位是km
# 注意传入lat和long的时候不要写反
# example: my_haversine(47.80844, -122.43159, 26.94168, -80.80621)
def my_haversine(lat1, long1, lat2, long2):
    lat1, long1, lat2, long2 = map(math.radians, [lat1, long1, lat2, long2])
    dlon = long2 - long1
    dlat = lat2 - lat1
    a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    c = 2 * math.asin(min(1, math.sqrt(a)))
    R = 6371.0088
    d = R * c
    return d

用google map的距离测算验证一组（粗略验证，因为google map不方便在测算距离的时候复制高精度经纬度）：

(47.80844, -122.43159)
到
(26.94168, -80.80621)

我的函数计算结果：

google map的距离测算：

再导入一个haversine包：🔗 [mapado/haversine: Calculate the distance between 2 points on Earth] https://github.com/mapado/haversine

跑10000组随机经纬度进行结果对比：

import math
import random
import haversine


# 注意返回的单位是km
# 注意传入lat和long的时候不要写反
# example: my_haversine(47.80844, -122.43159, 26.94168, -80.80621)
def my_haversine(lat1, long1, lat2, long2):
    lat1, long1, lat2, long2 = map(math.radians, [lat1, long1, lat2, long2])
    dlon = long2 - long1
    dlat = lat2 - lat1
    a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    c = 2 * math.asin(min(1, math.sqrt(a)))
    R = 6371.0088
    d = R * c
    return d


for i in range(0, 10 ** 4):
    rand_long1 = random.uniform(-180, 180)
    rand_long2 = random.uniform(-180, 180)
    rand_lat1 = random.uniform(-90, 90)
    rand_lat2 = random.uniform(-90, 90)
    my_haversine_dist = my_haversine(lat1=rand_lat1, long1=rand_long1, lat2=rand_lat2, long2=rand_long2)
    standard_haversine_dist = haversine.haversine((rand_lat1, rand_long1), (rand_lat2, rand_long2))
    if abs(my_haversine_dist - standard_haversine_dist) > 10 ** (-10):
        print(f'not equal: {abs(my_haversine_dist - standard_haversine_dist)}')
        exit(1)
print('pass')

通过了10000组随机测试，所以认为上面自己写的 my_haversine 没有什么问题

(long1, lat1), (long2, lat2)连线中任意一点的坐标计算

两种不同的解释方法

假设有一个参数用来表示两点间连线的任意一点，求这个点的经纬度坐标（用long和lat表示）

这个“连线”有2种解释：（下面假设这两个点在地球上距离相当远，已经不能忽略球体曲面的影响了）

解释1：在球体地图上表现为曲线，在展开地图上表现为直线

（可视化网站是geojson.io，geojson文件：解释1的geojson）

解释2：在球体地图上表现为直线，在展开地图上表现为曲线

（geojson文件：解释2的geojson）

我也不知道哪一种才是GIS业界“默认”的，只能把这两种解释的代码都贴出来。

解释1（by straight-line distance）

解释1：在球体地图上表现为曲线，在展开地图上表现为直线

实现代码和原理非常非常简单。就是把经纬度的坐标直接当成平面直角坐标系里的(x, y)，然后用线性代数的向量计算就可以搞定。

如果调用 shapely 包里的 line.interpolate 计算，也会使用这种“简单易懂“的2D平面坐标系进行计算。

shapely:

def middle_shapely(alpha, long1, lat1, long2, lat2):
    line = LineString([Point(long1, lat1), Point(long2, lat2)])
    # Calculate the midpoint with alpha = 0.5
    midpoint = line.interpolate(alpha, normalized=True)
    # Display the midpoint's coordinates
    return midpoint.x, midpoint.y  # x is longitude, y is latitude

自己写的：

def middle_my(alpha, long1, lat1, long2, lat2):
    return long1 + (long2 - long1) * alpha, lat1 + (lat2 - lat1) * alpha

对比代码：

from shapely.geometry import LineString, Point
import geojson
import numpy as np
import os

lat1, long1 = 47.85391930413978, -121.7528973290274
lat2, long2 = 26.382806911861756, -80.84341092679324


def middle_shapely(alpha, long1, lat1, long2, lat2):
    line = LineString([Point(long1, lat1), Point(long2, lat2)])
    # Calculate the midpoint with alpha = 0.5
    midpoint = line.interpolate(alpha, normalized=True)
    # Display the midpoint's coordinates
    return midpoint.x, midpoint.y  # x is longitude, y is latitude


def middle_my(alpha, long1, lat1, long2, lat2):
    return long1 + (long2 - long1) * alpha, lat1 + (lat2 - lat1) * alpha


# shapely
points = [(x, y) for x, y in [middle_shapely(alpha, long1, lat1, long2, lat2) for alpha in np.linspace(0, 1, 101)]]
features = [geojson.Feature(geometry=geojson.Point((x, y))) for x, y in points]
feature_collection = geojson.FeatureCollection(features)
geojson_data = geojson.dumps(feature_collection, indent=2)
with open("points_shapely.geojson", "w") as file:
    file.write(geojson_data)

# 自己写的
points = [(x, y) for x, y in [middle_my(alpha, long1, lat1, long2, lat2) for alpha in np.linspace(0, 1, 101)]]
features = [geojson.Feature(geometry=geojson.Point((x, y))) for x, y in points]
feature_collection = geojson.FeatureCollection(features)
geojson_data = geojson.dumps(feature_collection, indent=2)
with open("points_my.geojson", "w") as file:
    file.write(geojson_data)

# 对比结果
os.system('diff ./points_shapely.geojson ./points_my.geojson')

对比结果为：两个geojson文件完全一致。

解释2（by geodesic）