detect_utm_zone.py 文件源码

def get_bbox(shapefile):
    """ Gets the bounding box of a shapefile in EPSG 4326.
        If shapefile is not in WGS84, bounds are reprojected.
    """
    driver = ogr.GetDriverByName('ESRI Shapefile')
    # 0 means read, 1 means write
    data_source = driver.Open(shapefile, 0)
    # Check to see if shapefile is found.
    if data_source is None:
        failure('Could not open %s' % (shapefile))
    else:
        layer = data_source.GetLayer()
        shape_bbox = layer.GetExtent()

        spatialRef = layer.GetSpatialRef()
        target = osr.SpatialReference()
        target.ImportFromEPSG(4326)

        # this check for non-WGS84 projections gets some false positives, but that's okay
        if target.ExportToProj4() != spatialRef.ExportToProj4():
            transform = osr.CoordinateTransformation(spatialRef, target)
            point1 = ogr.Geometry(ogr.wkbPoint)
            point1.AddPoint(shape_bbox[0], shape_bbox[2])
            point2 = ogr.Geometry(ogr.wkbPoint)
            point2.AddPoint(shape_bbox[1], shape_bbox[3])
            point1.Transform(transform)
            point2.Transform(transform)
            shape_bbox = [point1.GetX(), point2.GetX(), point1.GetY(), point2.GetY()]

        return shape_bbox