shp平滑处理

　　在做图像数据处理时，经常会有栅格数据转矢量数据的操作，转换后的矢量文件会存在锯齿状边缘，不太美观，因此常常需要对矢量（shp）文件做平滑处理。

 

1 利用arcgis实现shp的平滑和简化

　　ArcToolbox / Cartography Tool / Generalization / Smooth Polygon，或者，制图工具 / 制图综合 / 平滑面。

 

2 基于GDAL矢量文件边界平滑

import sys
from osgeo import ogr, gdal, osr


#平滑函数
def smooth(shp,out,name,smooth_size):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    dataSource = driver.Open(shp, 0)  # 0是只读，1可写
    if dataSource is None:
        print('could not open')
        sys.exit(1)
    # 获取图层
    layer = dataSource.GetLayer(0)
    t = int(layer.GetFeatureCount())
    drv = ogr.GetDriverByName('ESRI Shapefile')
    Polygon = drv.CreateDataSource(out)
    #  oLayer = Polygon.CreateLayer(name, srs=prj, geom_type=ogr.wkbMultiPolygon)
    oLayer = Polygon.CreateLayer(name, layer.GetSpatialRef(), geom_type=ogr.wkbPolygon)
    oFieldID = ogr.FieldDefn("ID", ogr.OFTInteger)  # 创建一个叫ID的整型属性
    oLayer.CreateField(oFieldID, 1)
    feature = ogr.Feature(oLayer.GetLayerDefn())
    ID=0
    for i in range(0, t):
        feat = layer.GetFeature(i)
        geom = feat.GetGeometryRef()
        ID = ID+1
        buffer = geom.Buffer(smooth_size).Buffer(-smooth_size)
        feature.SetGeometry(buffer)  ## 设置面
        feature.SetField(0, ID)
        oLayer.CreateFeature(feature)


if __name__ == '__main__':
    #  shp路径
    shp=r"E:\成果\landsat8成果\LC08_L2SP_141039_20170214_20200905_02_T1.shp"
    #  输出路径
    out=r"E:\成果\修复"
    #  输出名称
    name=r"20170214"
    #  平滑尺度，投影为经纬度选择0.001为宜，投影为坐标选择100为宜。
    smooth_size = 0.001
    smooth(shp, out, name, smooth_size)