使用GDAL在Python脚本中简化shapefile的几何结构

2024-06-17 14:57:49 发布

您现在位置:Python中文网/ 问答频道 /正文

我正在尝试运行一个Python脚本,它可以简化shapefile中多边形的几何体。它似乎可以工作,但当我试图在像QGIS或gvSIG这样的GIS软件中渲染它时,shapefile不能显示几何图形。只有当文件大(超过60MB)时才会发生这种情况。在

我知道渲染它需要一些时间,但在某些时候软件必须渲染它。我不知道你对此有什么解释。在

下面是Python脚本:

import os
import time
import multiprocessing as mp
from osgeo import ogr
from osgeo import gdal

def multipoly2poly(in_lyr, out_lyr):
    for in_feat in in_lyr:
        geom = in_feat.GetGeometryRef()
        if geom.GetGeometryName() == 'MULTIPOLYGON':
            for geom_part in geom:
                x = geom_part.Simplify(0.05, preserve_topology=False)
                addPolygon(x.ExportToWkb(), out_lyr)
        else:
            addPolygon(geom.ExportToWkb(), out_lyr)

def addPolygon(simplePolygon, out_lyr):
    featureDefn = out_lyr.GetLayerDefn()
    polygon = ogr.CreateGeometryFromWkb(simplePolygon)
    out_feat = ogr.Feature(featureDefn)
    out_feat.SetGeometry(polygon)
    out_lyr.CreateFeature(out_feat)

def callback(result):
    print result

def processFolder(folderPath, driver):
    outputDir = folderPath+"/output"
    if not os.path.exists(outputDir):
        os.mkdir(outputDir)

    pool = mp.Pool(processes=4)
    shapefiles = [f for f in os.listdir(folderPath) if f.endswith(".shp")]
    time1 = time.time()
    for shapefile in shapefiles:
        processFile(shapefile, folderPath, outputDir, driver)

    time2 = time.time()
    print 'Total processing took %0.3f min' % ((time2-time1)/60.0)

def processFile(filePath, inputDir, outputDir, driver):
    print filePath    
    in_ds = driver.Open(inputDir+'/'+filePath, 0)
    in_lyr = in_ds.GetLayer()
    outputshp = outputDir+'/'+filePath
    print outputshp
    if os.path.exists(outputshp):
        print 'Shape %s exists, doing nothing' % (outputshp)
        #driver.DeleteDataSource(outputshp)
        return

    out_ds = driver.CreateDataSource(outputshp)
    out_lyr = out_ds.CreateLayer('poly', geom_type=ogr.wkbPolygon)

    time1 = time.time()
    multipoly2poly(in_lyr, out_lyr)
    time2 = time.time()
    print '%s took %0.3f min' % (filePath, (time2-time1)/60.0)
    return time2 - time1

gdal.UseExceptions()
driver = ogr.GetDriverByName('ESRI Shapefile')

processFolder(".", driver)

Tags: inimporttimeosdefdriveroutoutputdir