多边形向内收缩指定距离的实现

import math
import numpy as np
def calc_cos_angle(a,b):                #计算两个向量的夹角余弦
    inner_ab=a[0]*b[0]+a[1]*b[1]
    abs_a=math.sqrt(a[0]**2+a[1]**2)
    abs_b=math.sqrt(b[0]**2+b[1]**2)
    return inner_ab/(abs_a*abs_b)

def calc_point(p1,p2,p3,d):
    p1,p2,p3=np.array(p1),np.array(p2),np.array(p3)
    a=p1-p2
    b=p3-p2
    cos_angle=calc_cos_angle(a,b)                           #计算向量夹角余弦
    sin_half_angle=math.sqrt((1-cos_angle)/2)               #计算半角正弦值
    l=d/sin_half_angle                                      #计算要计算的点与中间点的距离
    abs_a=math.sqrt(a[0]**2+a[1]**2)                        #归一化向量a,b
    abs_b=math.sqrt(b[0]**2+b[1]**2)
    vector=a/abs_a+b/abs_b                                 #要求的点在角平分线上，即在单位化的向量a,b的和的方向上
    abs_vector=math.sqrt(vector[0]**2+vector[1]**2)        #归一化角平分线向量
    vector=vector/abs_vector*l                             #将角平分线向量的长度伸缩为l
    return p2+vector,p2-vector                             #计算两个可能的候选点


def shrink_polygon(polygon,d):
    x, y = polygon.exterior.coords.xy
    points=list(zip(x,y))
    points=points[0:-1]
    new_points=[]
    for i in range(len(points)):
        if i==0:
            p1=points[-1]
            p2=points[0]
            p3=points[1]
        elif i==len(points)-1:
            p1=points[-2]
            p2=points[-1]
            p3=points[0]
        else:
            p1=points[i-1]
            p2=points[i]
            p3=points[i+1]
        possible_point1,possible_point2=calc_point(p1,p2,p3,d)
        geo_possible_point1=geo.Point(possible_point1)
        geo_possible_point2=geo.Point(possible_point2)
        if geo_possible_point1.within(polygon):               #两个候选点在多边形内部时才是真正想要的点
            new_points.append(possible_point1)
        elif geo_possible_point2.within(polygon):
            new_points.append(possible_point2)
        else:                                                 #两点均落在多边形外部，则丢弃
            pass
    return geo.Polygon(new_points)

import shapely.geometry as geo

polygon=geo.Polygon([(0,0),(1.6,8),(1.25,2.25),(4,0),(1,-1)])
ax=geopandas.GeoSeries(polygon).plot()
new_polygon=shrink_polygon(polygon,0.1)
geopandas.GeoSeries(new_polygon).plot(ax=ax,color='r')
plt.show()

使用自带函数的实现：

#几何体缩放方法
def shrink1(polygon,distance,mode):
    if mode==0:
        polygon_bounds=geopandas.GeoSeries(polygon).boundary.iloc[0]
        polygon_bounds=polygon_bounds.parallel_offset(distance=distance, side='right', resolution=1, join_style=2, mitre_limit=5.0)
        return geo.Polygon(polygon_bounds.coords)
    else:
        return polygon.buffer(distance=-distance,resolution=1, cap_style=3, join_style=2, mitre_limit=5.0)

for i in range(len(gdf)):
    old=gdf.loc[i,'geometry']
    new=shrink1(old,15,1)
    ax=geopandas.GeoSeries(old).boundary.plot(color='b')
    geopandas.GeoSeries(new).boundary.plot(ax=ax,color='r')
    plt.show()