geolib.py 文件源码

def shp2geom(shp_fn):
    """Extract geometries from input shapefile

    Need to handle multi-part geom: http://osgeo-org.1560.x6.nabble.com/Multipart-to-singlepart-td3746767.html
    """
    ds = ogr.Open(shp_fn)
    lyr = ds.GetLayer()
    srs = lyr.GetSpatialRef()
    lyr.ResetReading()
    geom_list = []
    for feat in lyr:
        geom = feat.GetGeometryRef()
        geom.AssignSpatialReference(srs)
        #Duplicate the geometry, or segfault
        #See: http://trac.osgeo.org/gdal/wiki/PythonGotchas
        #g = ogr.CreateGeometryFromWkt(geom.ExportToWkt())
        #g.AssignSpatialReference(srs)
        g = geom_dup(geom)
        geom_list.append(g)
    #geom = ogr.ForceToPolygon(' '.join(geom_list))    
    #Dissolve should convert multipolygon to single polygon 
    #return geom_list[0]
    ds = None
    return geom_list