POI经纬度坐标常见用法整理（Python版）

本文章整理了经纬度坐标的相关用法，包括与geohash编码互相转换、火星坐标系转换算法、求两个坐标的距离、求围栏面积以及判断poi坐标是否在围栏里，后续再有经纬度相关的也会补充到这里~~芝士就是力量！power！！

目录

1.根据经纬度坐标获得geohash

2.GCJ-02 和 WGS 84 坐标系互相转换

GCJ-02 和 WGS 84 坐标系区别

火星坐标系转换算法

3.求两个坐标的距离

方法一

方法二

4.求围栏面积

5.判断poi坐标是否在围栏里

---

1.根据经纬度坐标获得geohash

输入不同的precision，可以获得不同长度的geohash，以geohash7为例：

import geohash2

# 将经度和纬度转换为 Geohash7 编码
def encode_geohash7(latitude, longitude):
    return geohash2.encode(latitude, longitude, precision=7)

逆操作：通过geohash编码转换为经纬度坐标

# 将 Geohash7 编码转换为经度和纬度
def decode_geohash7(geohash7):
    lat, lon = geohash2.decode(geohash7)
    return lat, lon

2.GCJ-02 和 WGS 84 坐标系互相转换

GCJ-02 和 WGS 84 坐标系区别

GCJ-02（国测局坐标系）和 WGS 84（World Geodetic System 1984）是两种不同的地理坐标系，用于表示地球上的位置。以下是它们之间的主要差别：

• 坐标偏移：
      WGS 84： 是一种全球通用的坐标系，由GPS使用。WGS 84 坐标通常被认为是地球上某一点的真实坐标。
      GCJ-02： 是中国国测局为保护国家地理信息安全而对 WGS 84 进行的一种加密处理。GCJ-02 对原始的 WGS 84 坐标进行了一定的偏移，使得在中国境内的地理位置不容易被外部解析。

• 使用范围：
      WGS 84： 广泛用于全球，是GPS设备和地理信息系统（GIS）等应用的标准坐标系。
      GCJ-02： 主要在中国境内使用，包括大部分在线地图服务如百度地图、高德地图等。

• 加密算法：
      WGS 84： 没有加密处理，坐标表示地球上的真实位置。
      GCJ-02： 使用了一种加密算法，通常称为"火星坐标系"。这个算法会对 WGS 84 坐标进行加密，以保护地理信息的安全。

• 坐标转换：
      由于 GCJ-02 的加密处理，WGS 84 坐标不能直接转换为 GCJ-02 坐标，反之亦然。因此，在进行坐标转换时，需要使用专门的算法，如火星坐标系转换算法。

火星坐标系转换算法

import math

#GCJ-02 坐标转换为 WGS 84 坐标
def gcj02_to_wgs84(gcj_lon, gcj_lat):
    a = 6378245.0
    ee = 0.00669342162296594323

    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 - ee * magic * magic
    sqrt_magic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrt_magic) * math.pi)
    dlon = (dlon * 180.0) / (a / sqrt_magic * math.cos(rad_lat) * math.pi)
    wgs_lat = gcj_lat - dlat
    wgs_lon = gcj_lon - dlon

    return wgs_lon, wgs_lat


#WGS 84 坐标转换为 GCJ-02 坐标
def wgs84_to_gcj02(wgs_lon, wgs_lat):
    a = 6378245.0
    ee = 0.00669342162296594323

    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 - ee * magic * magic
    sqrt_magic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrt_magic) * math.pi)
    dlon = (dlon * 180.0) / (a / sqrt_magic * math.cos(rad_lat) * math.pi)
    gcj_lat = wgs_lat + dlat
    gcj_lon = wgs_lon + dlon

    return gcj_lon, gcj_lat

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

3.求两个坐标的距离

介绍两种方法

方法一

from math import radians, sin, cos, sqrt, atan2

def haversine(lat1, lon1, lat2, lon2):
    # 将经纬度转换为弧度
    lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2])

    # Haversine 公式
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
    c = 2 * atan2(sqrt(a), sqrt(1 - a))

    # 地球半径（单位：千米）
    R = 6371.0

    # 计算距离
    distance = R * c

    return distance

# 示例：计算两个坐标点的距离
lat1, lon1 = 45.7890, 123.4560    
lat2, lon2 = 45.7899, 123.4560

distance_km = haversine(lat1, lon1, lat2, lon2)

# 打印结果
print(f"The distance between the two points is {distance_km*1000:.2f} m.")

方法二

from geopy.distance import geodesic

lat1, lon1 = (45.7890, 123.4560)  # 示例坐标1
lat2, lon2 = (45.7899, 123.4560)  # 示例坐标2


# 计算两个坐标之间的距离
distance_km = geodesic((lat1, lon1), (lat2, lon2)).kilometers

# 打印结果
print(f"The distance between the two points is {distance_km*1000:.2f} m.")

4.求围栏面积

尝试了多种求围栏面积的方法，对比后发现以下方法是比较准确的

参考来源：根据经纬度点求面积_经纬度 面积-CSDN博客

#计算围栏面积函数
def ComputeArea(temp):
    arr_len = len(temp)
    if arr_len < 3:
        return 0.0

    s = temp[0][1] * (temp[arr_len -1][0]-temp[1][0])
    for i in range(1,arr_len):
        s += temp[i][1] * (temp[i-1][0] - temp[(i+1)%arr_len][0])
    return round(math.fabs(s/2)*9101160000.085981,6)


# 围栏
polygon='113.986907,22.595623;113.988259,22.595658;113.988224,22.595045;113.986888,22.595036'

#处理围栏格式
polygon_gcj =[]
for lon_lat in polygon.split(';'):
    lon=float(lon_lat.split(',')[0])
    lat=float(lon_lat.split(',')[1])
    polygon_gcj +=[(lon,lat)]

print(polygon_gcj)

#计算面积
area=ComputeArea(polygon_gcj)

print("围栏面积是：",area)

5.判断poi坐标是否在围栏里

简单版

from shapely.geometry import Point, Polygon

# POI坐标
poi = (113.988188,22.595069)
polygon='112.964302,28.225389;112.974781,28.226108;112.970758,28.209185;112.966764,28.190226;112.964703,28.164677;112.96312,28.157887;112.960295,28.149132;112.958914,28.139693;112.953478,28.139737;112.957256,28.163486;112.959139,28.208019'

# 创建Shapely的Point和Polygon对象
poi = Point(poi)
fence = Polygon(polygon)

# 判断POI是否在围栏内
is_inside = poi.within(fence)

if is_inside:
    print("POI在围栏内")
else:
    print("POI不在围栏内")

数组版

from shapely.geometry import Point, Polygon
'''
假设有一个数组data，包含longitude，latitude，polyline列，
需要判断(longitude,latitude)是否在polyline里
'''

#单个封闭围栏格式处理，经纬度对用";"隔开
def lon_lat_list(fence_i):
    fence_l = []
    for lon_lat in fence_i.split(';'):
        lon = float(lon_lat.split(',')[0])
        lat = float(lon_lat.split(',')[1])
        fence_l += [(lon, lat)]
    return fence_l


#这里考虑了polyline包含多个封闭围栏的情况，每个封闭围栏间用"|"隔开
def polyline_to_fence(polyline):
    fence=[]
    if '|' in polyline:
        for fence_i in polyline.split('|'):
            fence_l=lon_lat_list(fence_i)
            fence.append(fence_l)
    else:
        fence_l = lon_lat_list(polyline)
        fence.append(fence_l)
    return fence

def is_in_fence(row):
    lon,lat = row['longitude'],row['latitude']
    poi=(lon,lat)

    # 创建Shapely的Point对象
    poi = Point(poi)
    
    #应用自定义函数处理围栏格式
    fence = polyline_to_fence(row['polyline_x'])

    count=0
    for fence_x in fence:
        # 创建Shapely的Polygon对象
        fence = Polygon(fence_x)

        # 判断POI是否在围栏内
        is_inside = poi.within(fence)
        
        #is_inside为布尔型：Ture or False
        if is_inside:
            count+=1

    #返回1表示在围栏，0表示不在围栏
    return count    


#apply应用数组的自定义函数
data['is_in_fence']=data.apply(is_in_fence,axis=1)

#输出在/不在围栏的数量
print(data['is_in_fence'].value_counts())