shputils.py 文件源码

def load_polygons_from_shapefile(filename, target_sr):
    """
    Loads the given shapefiles, reprojects the polygons it contains in the
    given target spatialreference.
    Returns:
        polygons: A list of polygons (list of points)
        attributes: A list of attributes (dict of strings)
    """
    shape = ogr.Open(filename)
    assert shape, "Couldn't open %s" % filename
    assert shape.GetLayerCount() == 1
    layer = shape.GetLayer()
    nfeatures = layer.GetFeatureCount()

    shape_sr = osr.SpatialReference()
    shape_sr.ImportFromWkt(layer.GetSpatialRef().ExportToWkt())

    # Transform from shape to image coordinates
    transform = osr.CoordinateTransformation(shape_sr, target_sr)

    # http://geoinformaticstutorial.blogspot.ch/2012/10/accessing-vertices-from-polygon-with.html
    polygons = []
    attributes = []
    for i in xrange(nfeatures):
        feature = layer.GetFeature(i)
        attr = feature.items()
        newattr = {}
        # some attributes contains unicode character. ASCIIfy everything
        # TODO: A bit brutal, but easy...
        for k, v in attr.items():
            newk = k.decode('ascii', errors='ignore')
            newattr[newk] = v
        attr = newattr

        geometry = feature.GetGeometryRef()
        assert geometry.GetGeometryName() == 'POLYGON'
        # A ring is a polygon in shapefiles
        ring = geometry.GetGeometryRef(0)
        assert ring.GetGeometryName() == 'LINEARRING'
        # The ring duplicates the last point, so for the polygon to be closed,
        # last point should equal first point
        npoints = ring.GetPointCount()
        points = [ring.GetPoint(i) for i in xrange(npoints)]
        points = [transform.TransformPoint(*p) for p in points]
        points = np.array(points)[:, :2]  # third column is elevation - discard
        # swap (lng, lat) to (lat, lng)
        points = points[:, ::-1]
        assert np.allclose(points[-1], points[0])
        polygons.append(points)
        attributes.append(attr)

    return polygons, attributes