shapely 经纬度相关计算

目录

判断经纬度是否落在轮廓里

获得轮廓的中心点

轮廓的中心点，及各个坐标串到中心点的最小距离

经纬度扩充

---

判断经纬度是否落在轮廓里

经纬度 lng，lat

轮廓 polygon，即经纬度坐标串，取值格式如下：

POLYGON ((113.757503 23.016089, 113.75749818472667 23.01599098285967, 113.7574837852804 23.015893909677985, 113.75745994033572 23.015798715322745, 113.7574268795325 23.015706316567634, 113.75738492126435 23.015617603263173))

is_in 得到的是 bool 值

from shapely import wkt

is_in = wkt.loads(polygon).contains(wkt.loads("POINT (" + str(lng) + ' ' + str(lat) + ')'))

获得轮廓的中心点

polygon 取值格式如第一个标题样例，在计算是必须转换为 shapely 能计算的对象

from shapely import wkt
from shapely.geometry import Point,polygon

polygon = wkt.loads(polygon) # 转换能计算的 shapely 对象
center = str(poly.centroid.x) + " " + str(poly.centroid.y)

或者

center = polygon.centroid

轮廓的中心点，及各个坐标串到中心点的最小距离

def def_projectpoly(project_poly):
    '''
    计算每条数据轮廓的中心点，及各个坐标串到中心点的最小距离
    '''
    # poly = project_poly['polygon'].iloc[0]
    def def_poly_min_distance(poly):
        poly_center = str(poly.centroid.x) + " " + str(poly.centroid.y)
        poly_wkt = poly.wkt.replace("POLYGON ","").replace("(","").replace(")","").replace("MULTI","")
        poly_wkt_list = poly_wkt.split(", ")
        # 计算坐标串各个组到中心点最短距离
        min_distance = 1000000
        for ind in poly_wkt_list:
            distance = geodistance(float(poly_center.split(" ")[0]),float(poly_center.split(" ")[1]),float(ind.split(" ")[0]),float(ind.split(" ")[1]))
            if distance < min_distance:
                min_distance = distance
        return min_distance

    # 将轮廓坐标串转换为 shape 对象
    project_poly['polygon'] = project_poly['polygon'].map(lambda x: wkt.loads(x))
    project_poly['min_distance'] = project_poly['polygon'].apply(lambda x:def_poly_min_distance(x))
    project_poly['poly_center'] = project_poly['polygon'].apply(lambda x:str(x.centroid.x) + " " + str(x.centroid.y))
    return project_poly


def geodistance(lng1, lat1, lng2, lat2):
    '''
    计算两两坐标的距离，在此是计算经纬度到轮廓中心点的距离
    '''
    # lng1,lat1,lng2,lat2 = (120.12802999999997,30.28708,115.86572000000001,28.7427)
    lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)])  # 经纬度转换成弧度
    dlon = lng2 - lng1
    dlat = lat2 - lat1
    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
    distance = 2 * asin(sqrt(a)) * 6371 * 1000  # 地球平均半径，6371km
    distance = round(distance / 1000, 3)
    return distance

经纬度扩充

扩充 50m，100m 写法如下

如果轮廓坐标串数量小于 10，则每个经纬度扩充 100m

from shapely import wkt
from shapely.geometry import Point,polygon

# 扩充范围在这里，0.0005 为 50m，0.001 为 100m
project['polygon'] = project['polygon'].apply(lambda x:x.buffer(0.0005) if len(str(x)) > 10 else '')## 增加在这里改

project['polygon'] = project.apply(lambda row:Point(row['LNG'],row['LAT']).buffer(0.001) if len(str(row['polygon'])) < 10 else row['polygon'],axis=1)