经纬度：距离计算、经纬度转utm、不同标准经纬度转换

1、根据经纬度计算距离

方案一：geodesic

geodesic方法返回的是两个经纬度坐标之间的最短地表距离，这是基于地球的实际椭球形状计算的。这也被称为大地线或测地线距离。
geodesic方法默认返回的距离单位是公里（km）。但是，你可以使用meters、miles、kilometers等属性来获取其他单位的距离。

from geopy.distance import geodesic
# 纬、经度坐标
coords_1 = (28.238832038269299, 112.99606908538851)
coords_2 = (28.240085702120774, 113.00430762487889)
# 计算距离
print(geodesic(coords_1, coords_2).meters)

方案二：distance

distance函数默认使用geodesic方法.

from geopy.distance import distance, great_circle
coords_1 = (28.238832038269299, 112.99606908538851)
coords_2 = (28.240085702120774, 113.00430762487889)
# distance函数默认使用geodesic方法
dist = distance(coords_1, coords_2).meters
print(dist)

# great_circle方法来计算两点之间的大圆距离。但是，由于地球不是完美的球体，great_circle方法的结果可能会稍有偏差。
dist2 = great_circle(coords_1, coords_2).meters
print(dist2)

方案三：haversine

"haversine"距离假设地球是一个完美的球体，然后使用haversine公式来计算两点之间的距离。这种方法简单易用，但是由于地球实际上是一个椭球体，所以"haversine"距离的结果可能不够精确。

# pip install haversine
from haversine import haversine
coords_1 = (28.238832038269299, 112.99606908538851)
coords_2 = (28.240085702120774, 113.00430762487889)
dist = haversine(coords_1, coords_2, unit='m')  # 默认单位为公里，如果需要米，可以使用 haversine(coords_1, coords_2, unit='m')
print(dist)

2、经纬度转UTM

def get_utm_zone(longitude):
    '''
     经纬度是在球面上的坐标，这意味着经纬度之间的距离并不是恒定的。
    例如，纬度每隔1度的距离在地球上几乎是恒定的（约111公里），但经度每隔1度的距离随着纬度的增加而减小。在赤道，经度每隔1度约为111公里，但在极地附近，经度每隔1度的距离接近于0。
    因此，如果你的数据集覆盖的地理范围很大，那么在经纬度上直接计算距离或面积可能会导致错误。为了避免这种情况，你可能需要先将经纬度转换为某种平面坐标系统（比如UTM），然后再进行计算。
    这个公式：
    1、用于把WGS84经纬度转换为UTM ZONE。
    2、使用的经度值应该是-180到180的范围，负值表示西经，正值表示东经。如果你的经度值是0到360的范围，你需要先转换为-180到180的范围。
    :param longitude: 经度
    :return:
    '''
    return int((longitude + 180) / 6) + 1

3、不同标准经纬度转换

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
'''
1、涉及到的坐标种类介绍：
（1）epsg:4326（WGS84经纬度）,是全球通用坐标系。
（2）中国国内使用的加密GCJ-02坐标。
（3）本项目使用的坐标：百度bd09ll坐标。地址描述转经纬度使用的api见：https://lbsyun.baidu.com/faq/api?title=webapi/guide/webservice-geocoding-base

2、坐标转换参考
参考1：https://blog.csdn.net/chenxu0136/article/details/119361975
参考2：https://blog.csdn.net/weixin_65167287/article/details/135241547
'''

import math

# 设置常量
r_pi = math.pi * 3000.0 / 180.0
la = 6378245.0  # 长半轴
ob = 0.00669342162296594323  # 扁率

def out_of_china(lon, lat):
    return lon < 72.004 or lon > 137.8347 or lat < 0.8293 or lat > 55.8271


def wgs84_to_gcj02(wgs_lon, wgs_lat):
    '''
    # WGS 84 坐标转换为 GCJ-02 坐标
    :param wgs_lon: WGS 84经度
    :param wgs_lat: WGS 84维度
    :return: 转换后的经纬度坐标，[经度, 纬度]
    '''

    def transform_lon(x, y):
        ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
        return ret

    def transform_lat(x, y):
        ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320.0 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
        return ret

    if out_of_china(wgs_lon, wgs_lat):
        return wgs_lon, wgs_lat

    dlat = transform_lat(wgs_lon - 105.0, wgs_lat - 35.0)
    dlon = transform_lon(wgs_lon - 105.0, wgs_lat - 35.0)
    rad_lat = wgs_lat / 180.0 * math.pi
    magic = math.sin(rad_lat)
    magic = 1 - ob * magic * magic
    sqrt_magic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((la * (1 - ob)) / (magic * sqrt_magic) * math.pi)
    dlon = (dlon * 180.0) / (la / sqrt_magic * math.cos(rad_lat) * math.pi)
    gcj_lat = wgs_lat + dlat
    gcj_lon = wgs_lon + dlon

    return gcj_lon, gcj_lat


def bd09_to_gcj02(lon_bd09, lat_bd09):
    '''
     bd09 -> gcj02
    :param lat_bd09: bd09的纬度
    :param lon_bd09: bd09的经度
    :return: 转换后坐标的列表形式 ，[经度, 纬度]
    '''
    m = lon_bd09 - 0.0065
    n = lat_bd09 - 0.006
    c = math.sqrt(m * m + n * n) - 0.00002 * math.sin(n * r_pi)
    theta = math.atan2(n, m) - 0.000003 * math.cos(m * r_pi)
    lon_gcj02 = c * math.cos(theta)
    lat_gcj02 = c * math.sin(theta)
    return lon_gcj02, lat_gcj02


def gcj02_to_wgs84(gcj_lon, gcj_lat):
    '''
    GCJ-02 坐标转换为 WGS 84 坐标
    :param gcj_lon: gcj02的经度
    :param gcj_lat: gcj02的纬度
    :return: 转换后坐标的列表形式，[经度, 纬度]
    '''

    def transform_lon(x, y):
        ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
        return ret

    def transform_lat(x, y):
        ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320.0 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
        return ret

    if out_of_china(gcj_lon, gcj_lat):
        return gcj_lon, gcj_lat

    dlat = transform_lat(gcj_lon - 105.0, gcj_lat - 35.0)
    dlon = transform_lon(gcj_lon - 105.0, gcj_lat - 35.0)
    rad_lat = gcj_lat / 180.0 * math.pi
    magic = math.sin(rad_lat)
    magic = 1 - ob * magic * magic
    sqrt_magic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((la * (1 - ob)) / (magic * sqrt_magic) * math.pi)
    dlon = (dlon * 180.0) / (la / sqrt_magic * math.cos(rad_lat) * math.pi)
    wgs_lat = gcj_lat - dlat
    wgs_lon = gcj_lon - dlon

    return wgs_lon, wgs_lat


def bd09_to_wgs84(lon_bd09, lat_bd09):
    '''
     百度坐标系BD09 转 国际标准坐标系 WGS84
    :param lon_bd09: bd09经度
    :param lat_bd09: bd09维度
    :return: 转换后坐标的列表形式，[经度, 纬度]
    '''
    # 先把百度坐标系(bd09)的经纬度转换为火星坐标系(中国国内的标准坐标系GCJ02)
    gcj02_lon, gcj02_lat = bd09_to_gcj02(lon_bd09, lat_bd09)
    # 然后把火星坐标系的坐标转换为国际通用坐标系WGS84
    return gcj02_to_wgs84(gcj02_lon, gcj02_lat)