Python osgeo.ogr 模块，wkbPolygon() 实例源码

def exporttogeojson(geojsonfilename, buildinglist):
    #
    # geojsonname should end with .geojson
    # building list should be list of dictionaries
    # list of Dictionaries {'ImageId': image_id, 'BuildingId': building_id, 'polyPix': poly,
    #                       'polyGeo': poly}
    # image_id is a string,
    # BuildingId is an integer,
    # poly is a ogr.Geometry Polygon
    #
    # returns geojsonfilename

    driver = ogr.GetDriverByName('geojson')
    if os.path.exists(geojsonfilename):
        driver.DeleteDataSource(geojsonfilename)
    datasource = driver.CreateDataSource(geojsonfilename)
    layer = datasource.CreateLayer('buildings', geom_type=ogr.wkbPolygon)
    field_name = ogr.FieldDefn("ImageId", ogr.OFTString)
    field_name.SetWidth(75)
    layer.CreateField(field_name)
    layer.CreateField(ogr.FieldDefn("BuildingId", ogr.OFTInteger))

    # loop through buildings
    for building in buildinglist:
        # create feature
        feature = ogr.Feature(layer.GetLayerDefn())
        feature.SetField("ImageId", building['ImageId'])
        feature.SetField("BuildingId", building['BuildingId'])
        feature.SetGeometry(building['polyPix'])

        # Create the feature in the layer (geojson)
        layer.CreateFeature(feature)
        # Destroy the feature to free resources
        feature.Destroy()

    datasource.Destroy()

    return geojsonfilename

def SetZ (poGeom, dfZ ):

    if poGeom is None:
        return

    eGType = wkbFlatten(poGeom.GetGeometryType())
    if eGType == ogr.wkbPoint:
        poGeom.SetPoint(0, poGeom.GetX(), poGeom.GetY(), dfZ)

    elif eGType == ogr.wkbLineString or \
         eGType == ogr.wkbLinearRing:
        for i in range(poGeom.GetPointCount()):
            poGeom.SetPoint(i, poGeom.GetX(i), poGeom.GetY(i), dfZ)

    elif eGType == ogr.wkbPolygon or \
         eGType == ogr.wkbMultiPoint or \
         eGType == ogr.wkbMultiLineString or \
         eGType == ogr.wkbMultiPolygon or \
         eGType == ogr.wkbGeometryCollection:
        for i in range(poGeom.GetGeometryCount()):
            SetZ(poGeom.GetGeometryRef(i), dfZ)

#**********************************************************************
#                         SetupTargetLayer()
#**********************************************************************

def SetZ (poGeom, dfZ ):

    if poGeom is None:
        return

    eGType = wkbFlatten(poGeom.GetGeometryType())
    if eGType == ogr.wkbPoint:
        poGeom.SetPoint(0, poGeom.GetX(), poGeom.GetY(), dfZ)

    elif eGType == ogr.wkbLineString or \
         eGType == ogr.wkbLinearRing:
        for i in range(poGeom.GetPointCount()):
            poGeom.SetPoint(i, poGeom.GetX(i), poGeom.GetY(i), dfZ)

    elif eGType == ogr.wkbPolygon or \
         eGType == ogr.wkbMultiPoint or \
         eGType == ogr.wkbMultiLineString or \
         eGType == ogr.wkbMultiPolygon or \
         eGType == ogr.wkbGeometryCollection:
        for i in range(poGeom.GetGeometryCount()):
            SetZ(poGeom.GetGeometryRef(i), dfZ)

#/************************************************************************/
#/*                         SetupTargetLayer()                           */
#/************************************************************************/

def buildTindex(rasterFolder, rasterExtention='.tif'):
    rasterList = glob.glob(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
    print(rasterList)

    print(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))

    memDriver = ogr.GetDriverByName('MEMORY')
    gTindex = memDriver.CreateDataSource('gTindex')
    srcImage = gdal.Open(rasterList[0])
    spat_ref = osr.SpatialReference()
    spat_ref.SetProjection(srcImage.GetProjection())
    gTindexLayer = gTindex.CreateLayer("gtindexlayer", spat_ref, geom_type=ogr.wkbPolygon)

    # Add an ID field
    idField = ogr.FieldDefn("location", ogr.OFTString)
    gTindexLayer.CreateField(idField)

    # Create the feature and set values
    featureDefn = gTindexLayer.GetLayerDefn()



    for rasterFile in rasterList:
        srcImage = gdal.Open(rasterFile)

        geoTrans, polyToCut, ulX, ulY, lrX, lrY = gT.getRasterExtent(srcImage)

        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(polyToCut)
        feature.SetField("location", rasterFile)
        gTindexLayer.CreateFeature(feature)
        feature = None


    return gTindex, gTindexLayer

def returnBoundBox(xOff, yOff, pixDim, inputRaster, targetSR='', pixelSpace=False):
    # Returns Polygon for a specific square defined by a center Pixel and
    # number of pixels in each dimension.
    # Leave targetSR as empty string '' or specify it as a osr.SpatialReference()
    # targetSR = osr.SpatialReference()
    # targetSR.ImportFromEPSG(4326)
    if targetSR == '':
        targetSR = osr.SpatialReference()
        targetSR.ImportFromEPSG(4326)
    xCord = [xOff - pixDim / 2, xOff - pixDim / 2, xOff + pixDim / 2,
             xOff + pixDim / 2, xOff - pixDim / 2]
    yCord = [yOff - pixDim / 2, yOff + pixDim / 2, yOff + pixDim / 2,
             yOff - pixDim / 2, yOff - pixDim / 2]

    ring = ogr.Geometry(ogr.wkbLinearRing)
    for idx in xrange(len(xCord)):
        if pixelSpace == False:
            geom = pixelToGeoCoord(xCord[idx], yCord[idx], inputRaster)
            ring.AddPoint(geom[0], geom[1], 0)
        else:
            ring.AddPoint(xCord[idx], yCord[idx], 0)

    poly = ogr.Geometry(ogr.wkbPolygon)
    if pixelSpace == False:
        poly.AssignSpatialReference(targetSR)

    poly.AddGeometry(ring)

    return poly

def getRasterExtent(srcImage):
    geoTrans = srcImage.GetGeoTransform()
    ulX = geoTrans[0]
    ulY = geoTrans[3]
    xDist = geoTrans[1]
    yDist = geoTrans[5]
    rtnX = geoTrans[2]
    rtnY = geoTrans[4]

    cols = srcImage.RasterXSize
    rows = srcImage.RasterYSize

    lrX = ulX + xDist * cols
    lrY = ulY + yDist * rows

    # Create ring
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(lrX, lrY)
    ring.AddPoint(lrX, ulY)
    ring.AddPoint(ulX, ulY)
    ring.AddPoint(ulX, lrY)
    ring.AddPoint(lrX, lrY)

    # Create polygon
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)

    return geoTrans, poly, ulX, ulY, lrX, lrY

def createPolygonFromCorners(lrX,lrY,ulX, ulY):
    # Create ring
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(lrX, lrY)
    ring.AddPoint(lrX, ulY)
    ring.AddPoint(ulX, ulY)
    ring.AddPoint(ulX, lrY)
    ring.AddPoint(lrX, lrY)

    # Create polygon
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)

    return poly

def evaluateLineStringPlane(geom, label='Airplane'):
    ring = ogr.Geometry(ogr.wkbLinearRing)

    for i in range(0, geom.GetPointCount()):
        # GetPoint returns a tuple not a Geometry
        pt = geom.GetPoint(i)
        ring.AddPoint(pt[0], pt[1])
    pt = geom.GetPoint(0)
    ring.AddPoint(pt[0], pt[1])
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)

    transform_WGS84_To_UTM, transform_UTM_To_WGS84, utm_cs = gT.createUTMTransform(geom)
    geom.Transform(transform_WGS84_To_UTM)
    pt0 = geom.GetPoint(0) # Tail
    pt1 = geom.GetPoint(1) # Wing
    pt2 = geom.GetPoint(2) # Nose
    pt3 = geom.GetPoint(3) # Wing
    Length = math.sqrt((pt2[0]-pt0[0])**2 + (pt2[1]-pt0[1])**2)
    Width = math.sqrt((pt3[0] - pt1[0])**2 + (pt3[1] - pt1[1])**2)
    Aspect = Length/Width
    Direction = (math.atan2(pt2[0]-pt0[0], pt2[1]-pt0[1])*180/math.pi) % 360


    geom.Transform(transform_UTM_To_WGS84)

    return [poly, Length, Width, Aspect, Direction]

def _get_vector_attributes(self,ds,attrs):
        logging.debug("[FileMan][_get_vector_attributes]")

        layer = ds.GetLayer()

        # extent
        extent = layer.GetExtent()
        logging.debug("[FileMan][_get_vector_attributes] Extent: %s" % repr(extent) )
        attrs["extent"] = (extent[0],extent[2],extent[1],extent[3])

        # features count
        attrs["features_count"] = layer.GetFeatureCount()

        # geom type
        ftype = layer.GetGeomType()
        if ftype == ogr.wkbPolygon:    #3
            attrs["type"] = "polygon"
        elif ftype == ogr.wkbPoint:    #1
            attrs["type"] = "point"
        elif ftype == ogr.wkbLineString:   #2
            attrs["type"] = "line"
        elif ftype == ogr.wkbPolygon25D:
            attrs["type"] = "polygon"
        else: 
            attrs["type"] = "none/unknown"  # 0 or anything else

        # srs
        sr = layer.GetSpatialRef()
        logging.debug("[FileMan][_get_vector_attributes] Spatial Reference from layer.GetSpatialRef() : %s" % repr(sr) )
        if sr:
            attrs["prj"]= self._get_prj(sr)
        else:
            attrs["prj"] = "unknown"

        # Done
        return attrs

def _check_src(self, src):
        """ Basic check of source elements (sequence of points or polygons).

            - array cast of source elements
            - create ogr_src datasource/layer holding src points/polygons
            - transforming source grid points/polygons to ogr.geometries
              on ogr.layer
        """
        ogr_src = io.gdal_create_dataset('Memory', 'out',
                                         gdal_type=gdal.OF_VECTOR)

        try:
            # is it ESRI Shapefile?
            ds_in, tmp_lyr = io.open_shape(src, driver=ogr.
                                           GetDriverByName('ESRI Shapefile'))
            ogr_src_lyr = ogr_src.CopyLayer(tmp_lyr, self._name)
            if self._srs is None:
                self._srs = ogr_src_lyr.GetSpatialRef()
        except IOError:
            # no ESRI shape file
            raise
        # all failed? then it should be sequence or numpy array
        except RuntimeError:
            src = np.array(src)
            # create memory datasource, layer and create features
            if src.ndim == 2:
                geom_type = ogr.wkbPoint
            # no Polygons, just Points
            else:
                geom_type = ogr.wkbPolygon
            fields = [('index', ogr.OFTInteger)]
            georef.ogr_create_layer(ogr_src, self._name, srs=self._srs,
                                    geom_type=geom_type, fields=fields)
            georef.ogr_add_feature(ogr_src, src, name=self._name)

        return ogr_src

def ogr_create_layer(ds, name, srs=None, geom_type=None, fields=None):
    """Creates OGR.Layer objects in gdal.Dataset object.

    .. versionadded:: 0.7.0

    Creates one OGR.Layer with given name in given gdal.Dataset object
    using given OGR.GeometryType and FieldDefinitions

    Parameters
    ----------
    ds : gdal.Dataset
        object
    name : string
        OGRLayer name
    srs : OSR.SpatialReference
        object
    geom_type : OGR GeometryType
        (eg. ogr.wkbPolygon)
    fields : list of 2 element tuples
        (strings, OGR.DataType) field name, field type

    Returns
    -------
    out : OGR.Layer
        object
    """
    if geom_type is None:
        raise TypeError("geometry_type needed")

    lyr = ds.CreateLayer(name, srs=srs, geom_type=geom_type)
    if fields is not None:
        for fname, fvalue in fields:
            lyr.CreateField(ogr.FieldDefn(fname, fvalue))

    return lyr

def bbox(coordinates, crs, outname=None, format='ESRI Shapefile', overwrite=True):
    """
    create a bounding box vector object or shapefile from coordinates and coordinate reference system
    coordinates must be provided in a dictionary containing numerical variables with names 'xmin', 'xmax', 'ymin' and 'ymax'
    the coordinate reference system can be in either WKT, EPSG or PROJ4 format
    """
    srs = crsConvert(crs, 'osr')

    ring = ogr.Geometry(ogr.wkbLinearRing)

    ring.AddPoint(coordinates['xmin'], coordinates['ymin'])
    ring.AddPoint(coordinates['xmin'], coordinates['ymax'])
    ring.AddPoint(coordinates['xmax'], coordinates['ymax'])
    ring.AddPoint(coordinates['xmax'], coordinates['ymin'])
    ring.CloseRings()

    geom = ogr.Geometry(ogr.wkbPolygon)
    geom.AddGeometry(ring)

    geom.FlattenTo2D()

    bbox = Vector(driver='Memory')
    bbox.addlayer('bbox', srs, ogr.wkbPolygon)
    bbox.addfield('id', width=4)
    bbox.addfeature(geom, 'id', 1)
    geom.Destroy()
    if outname is None:
        return bbox
    else:
        bbox.write(outname, format, overwrite)

def export_shp_triangles(self, outshp):
        from osgeo import ogr

        # Now convert it to a shapefile with OGR
        driver = ogr.GetDriverByName('Esri Shapefile')
        ds = driver.CreateDataSource(outshp)
        layer = ds.CreateLayer('', None, ogr.wkbPolygon)
        # Add one attribute
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
        defn = layer.GetLayerDefn()

        ## If there are multiple geometries, put the "for" loop here
        for id, poly in self.triangles.iteritems():
            # Create a new feature (attribute and geometry)
            feat = ogr.Feature(defn)
            feat.SetField('id', int(id))

            # Make a geometry, from Shapely object
            geom = ogr.CreateGeometryFromWkb(poly.wkb)
            feat.SetGeometry(geom)

            layer.CreateFeature(feat)
            feat = geom = None  # destroy these

        # Save and close everything
        ds = layer = feat = geom = None

def export_shp_hull(self, outshp):
        from osgeo import ogr

        # Now convert it to a shapefile with OGR
        driver = ogr.GetDriverByName('Esri Shapefile')
        ds = driver.CreateDataSource(outshp)
        layer = ds.CreateLayer('', None, ogr.wkbPolygon)
        # Add one attribute
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
        defn = layer.GetLayerDefn()

        ## If there are multiple geometries, put the "for" loop here
        for id, poly in self.hull_polygons.iteritems():
            # Create a new feature (attribute and geometry)
            feat = ogr.Feature(defn)
            feat.SetField('id', int(id))

            # Make a geometry, from Shapely object
            geom = ogr.CreateGeometryFromWkb(poly.wkb)
            feat.SetGeometry(geom)

            layer.CreateFeature(feat)
            feat = geom = None  # destroy these

        # Save and close everything
        ds = layer = feat = geom = None

def get_convex_hull(shape_name, destination_directory):
    '''
    This method will read all objects from a shape file and create a new one
    with the convex hull of all the geometry points of the first.
    '''
    driver = ogr.GetDriverByName(str('ESRI Shapefile'))
    shape = driver.Open(shape_name, 0)
    layer = shape.GetLayer()
    layer_name = layer.GetName()
    spatial_reference = layer.GetSpatialRef()
    prefix = get_basename(shape_name)
    output_name = create_filename(destination_directory, '%s-hull.shp' % prefix)
    geometries = ogr.Geometry(ogr.wkbGeometryCollection)
    for feature in layer:
        geometries.AddGeometry(feature.GetGeometryRef())
    if os.path.exists(output_name):
        driver.DeleteDataSource(output_name)
    datasource = driver.CreateDataSource(output_name)
    out_layer = datasource.CreateLayer(str('states_convexhull'), spatial_reference, geom_type=ogr.wkbPolygon)
    out_layer.CreateField(ogr.FieldDefn(str('id'), ogr.OFTInteger))
    featureDefn = out_layer.GetLayerDefn()
    feature = ogr.Feature(featureDefn)
    feature.SetGeometry(geometries.ConvexHull())
    feature.SetField(str('id'), 1)
    out_layer.CreateFeature(feature)
    shape.Destroy()
    datasource.Destroy()

def createBoxFromLine(tmpGeom, ratio=1, halfWidth=-999, transformRequired=True, transform_WGS84_To_UTM='', transform_UTM_To_WGS84=''):
    # create Polygon Box Oriented with the line

    if transformRequired:
        if transform_WGS84_To_UTM == '':
            transform_WGS84_To_UTM, transform_UTM_To_WGS84 = createUTMTransform(tmpGeom)

        tmpGeom.Transform(transform_WGS84_To_UTM)


    # calculatuate Centroid
    centroidX, centroidY, centroidZ = tmpGeom.Centroid().GetPoint()
    lengthM = tmpGeom.Length()
    if halfWidth ==-999:
        halfWidth = lengthM/(2*ratio)

    envelope=tmpGeom.GetPoints()
    cX1 = envelope[0][0]
    cY1 = envelope[0][1]
    cX2 = envelope[1][0]
    cY2 = envelope[1][1]
    angRad = math.atan2(cY2-cY1,cX2-cX1)

    d_X = math.cos(angRad-math.pi/2)*halfWidth
    d_Y = math.sin(angRad-math.pi/2)*halfWidth

    e_X = math.cos(angRad+math.pi/2)*halfWidth
    e_Y = math.sin(angRad+math.pi/2)*halfWidth

    ring = ogr.Geometry(ogr.wkbLinearRing)

    ring.AddPoint(cX1+d_X, cY1+d_Y)
    ring.AddPoint(cX1+e_X, cY1+e_Y)
    ring.AddPoint(cX2+e_X, cY2+e_Y)
    ring.AddPoint(cX2+d_X, cY2+d_Y)
    ring.AddPoint(cX1+d_X, cY1+d_Y)
    polyGeom = ogr.Geometry(ogr.wkbPolygon)
    polyGeom.AddGeometry(ring)
    areaM = polyGeom.GetArea()

    if transformRequired:
        tmpGeom.Transform(transform_UTM_To_WGS84)
        polyGeom.Transform(transform_UTM_To_WGS84)


    return (polyGeom, areaM, angRad, lengthM)

def geoPolygonToPixelPolygonWKT(geom, inputRaster, targetSR, geomTransform, breakMultiPolygonGeo=True,
                                pixPrecision=2):
    # Returns Pixel Coordinate List and GeoCoordinateList

    polygonPixBufferList = []
    polygonPixBufferWKTList = []
    polygonGeoWKTList = []
    if geom.GetGeometryName() == 'POLYGON':
        polygonPix = ogr.Geometry(ogr.wkbPolygon)
        for ring in geom:
            # GetPoint returns a tuple not a Geometry
            ringPix = ogr.Geometry(ogr.wkbLinearRing)

            for pIdx in xrange(ring.GetPointCount()):
                lon, lat, z = ring.GetPoint(pIdx)
                xPix, yPix = latlon2pixel(lat, lon, inputRaster, targetSR, geomTransform)

                xPix = round(xPix, pixPrecision)
                yPix = round(yPix, pixPrecision)
                ringPix.AddPoint(xPix, yPix)

            polygonPix.AddGeometry(ringPix)
        polygonPixBuffer = polygonPix.Buffer(0.0)
        polygonPixBufferList.append([polygonPixBuffer, geom])

    elif geom.GetGeometryName() == 'MULTIPOLYGON':

        for poly in geom:
            polygonPix = ogr.Geometry(ogr.wkbPolygon)
            for ring in poly:
                # GetPoint returns a tuple not a Geometry
                ringPix = ogr.Geometry(ogr.wkbLinearRing)

                for pIdx in xrange(ring.GetPointCount()):
                    lon, lat, z = ring.GetPoint(pIdx)
                    xPix, yPix = latlon2pixel(lat, lon, inputRaster, targetSR, geomTransform)

                    xPix = round(xPix, pixPrecision)
                    yPix = round(yPix, pixPrecision)
                    ringPix.AddPoint(xPix, yPix)

                polygonPix.AddGeometry(ringPix)
            polygonPixBuffer = polygonPix.Buffer(0.0)
            if breakMultiPolygonGeo:
                polygonPixBufferList.append([polygonPixBuffer, poly])
            else:
                polygonPixBufferList.append([polygonPixBuffer, geom])

    for polygonTest in polygonPixBufferList:
        if polygonTest[0].GetGeometryName() == 'POLYGON':
            polygonPixBufferWKTList.append([polygonTest[0].ExportToWkt(), polygonTest[1].ExportToWkt()])
        elif polygonTest[0].GetGeometryName() == 'MULTIPOLYGON':
            for polygonTest2 in polygonTest[0]:
                polygonPixBufferWKTList.append([polygonTest2.ExportToWkt(), polygonTest[1].ExportToWkt()])

    return polygonPixBufferWKTList

def convertLabelStringToPoly(shapeFileSrc, outGeoJSon, labelType='Airplane'):

        shapeSrc = ogr.Open(shapeFileSrc)
        source_layer = shapeSrc.GetLayer()
        source_srs = source_layer.GetSpatialRef()
        # Create the output Layer
        outDriver = ogr.GetDriverByName("geojson")
        if os.path.exists(outGeoJSon):
            outDriver.DeleteDataSource(outGeoJSon)


        outDataSource = outDriver.CreateDataSource(outGeoJSon)
        outLayer = outDataSource.CreateLayer("groundTruth", source_srs, geom_type=ogr.wkbPolygon)
        # Add input Layer Fields to the output Layer
        inLayerDefn = source_layer.GetLayerDefn()
        for i in range(0, inLayerDefn.GetFieldCount()):
            fieldDefn = inLayerDefn.GetFieldDefn(i)
            outLayer.CreateField(fieldDefn)
        outLayer.CreateField(ogr.FieldDefn("Length_m", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("Width_m", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("Aspect(L/W)", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("compassDeg", ogr.OFTReal))

        outLayerDefn = outLayer.GetLayerDefn()
        for inFeature in source_layer:

            outFeature = ogr.Feature(outLayerDefn)

            for i in range(0, inLayerDefn.GetFieldCount()):
                outFeature.SetField(inLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))

            geom = inFeature.GetGeometryRef()
            if labelType == 'Airplane':
                poly, Length, Width, Aspect, Direction = evaluateLineStringPlane(geom, label='Airplane')
            elif labelType == 'Boat':
                poly, Length, Width, Aspect, Direction = evaluateLineStringBoat(geom, label='Boat')

            outFeature.SetGeometry(poly)
            outFeature.SetField("Length_m", Length)
            outFeature.SetField("Width_m", Width)
            outFeature.SetField("Aspect(L/W)", Aspect)
            outFeature.SetField("compassDeg", Direction)

            outLayer.CreateFeature(outFeature)

def polygon_to_shapefile(polygons, poly_sr, fname, fields_defs=None,
                         poly_attrs=None):
    """
    Write a set of polygons to a shapefile
    Args:
        polygons: a list of (lat, lng) tuples
        fields_defs: The list of fields for those polygons, as a tuple
                     (name, ogr type) for each field. For example :
                     [('field1', ogr.OFTReal), ('field2', ogr.OFTInteger)]
        poly_attrs: A list of dict containing the attributes for each polygon
                        [{'field1' : 1.0, 'field2': 42},
                         {'field1' : 3.0, 'field2': 60}]
    """
    shp_driver = ogr.GetDriverByName("ESRI Shapefile")
    out_ds = shp_driver.CreateDataSource(fname)
    assert out_ds is not None, "Failed to create temporary %s" % fname
    out_layer = out_ds.CreateLayer(fname, poly_sr, geom_type=ogr.wkbPolygon)

    has_attrs = fields_defs is not None
    if has_attrs:
        attrs_name = []
        for field_name, field_type in fields_defs:
            out_layer.CreateField(ogr.FieldDefn(field_name, field_type))
            attrs_name.append(field_name)

    layer_defn = out_layer.GetLayerDefn()
    for i, poly in enumerate(polygons):
        ring = ogr.Geometry(ogr.wkbLinearRing)
        # gdal uses the (x, y) convention => (lng, lat)
        for point in poly:
            ring.AddPoint(point[1], point[0])
        ring.AddPoint(poly[0][1], poly[0][0])  # re-add the start to close
        p = ogr.Geometry(ogr.wkbPolygon)
        p.AddGeometry(ring)

        out_feature = ogr.Feature(layer_defn)
        out_feature.SetGeometry(p)

        if has_attrs:
            attrs = poly_attrs[i]
            for field_name in attrs_name:
                out_feature.SetField(field_name, attrs[field_name])

        out_layer.CreateFeature(out_feature)

    out_feature.Destroy()
    out_ds.Destroy()

def DumpReadableGeometry( poGeometry, pszPrefix, options ):

    if pszPrefix == None:
        pszPrefix = ""

    if 'DISPLAY_GEOMETRY' in options and EQUAL(options['DISPLAY_GEOMETRY'], 'SUMMARY'):

        line = ("%s%s : " % (pszPrefix, poGeometry.GetGeometryName() ))
        eType = poGeometry.GetGeometryType()
        if eType == ogr.wkbLineString or eType == ogr.wkbLineString25D:
            line = line + ("%d points" % poGeometry.GetPointCount())
            print(line)
        elif eType == ogr.wkbPolygon or eType == ogr.wkbPolygon25D:
            nRings = poGeometry.GetGeometryCount()
            if nRings == 0:
                line = line + "empty"
            else:
                poRing = poGeometry.GetGeometryRef(0)
                line = line + ("%d points" % poRing.GetPointCount())
                if nRings > 1:
                    line = line + (", %d inner rings (" % (nRings - 1))
                    for ir in range(0,nRings-1):
                        if ir > 0:
                            line = line + ", "
                        poRing = poGeometry.GetGeometryRef(ir+1)
                        line = line + ("%d points" % poRing.GetPointCount())
                    line = line + ")"
            print(line)

        elif eType == ogr.wkbMultiPoint or \
            eType == ogr.wkbMultiPoint25D or \
            eType == ogr.wkbMultiLineString or \
            eType == ogr.wkbMultiLineString25D or \
            eType == ogr.wkbMultiPolygon or \
            eType == ogr.wkbMultiPolygon25D or \
            eType == ogr.wkbGeometryCollection or \
            eType == ogr.wkbGeometryCollection25D:

                line = line + "%d geometries:" % poGeometry.GetGeometryCount()
                print(line)
                for ig in range(poGeometry.GetGeometryCount()):
                    subgeom = poGeometry.GetGeometryRef(ig)
                    from sys import version_info
                    if version_info >= (3,0,0):
                        exec('print("", end=" ")')
                    else:
                        exec('print "", ')
                    DumpReadableGeometry( subgeom, pszPrefix, options)
        else:
            print(line)

    elif 'DISPLAY_GEOMETRY' not in options or EQUAL(options['DISPLAY_GEOMETRY'], 'yes') \
            or EQUAL(options['DISPLAY_GEOMETRY'], 'WKT'):

        print("%s%s" % (pszPrefix, poGeometry.ExportToWkt() ))

    return

def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere):
    poGeom = None

    poDS = ogr.Open( pszDS, False )
    if poDS is None:
        return None

    if pszSQL is not None:
        poLyr = poDS.ExecuteSQL( pszSQL, None, None )
    elif pszLyr is not None:
        poLyr = poDS.GetLayerByName(pszLyr)
    else:
        poLyr = poDS.GetLayer(0)

    if poLyr is None:
        print("Failed to identify source layer from datasource.")
        poDS.Destroy()
        return None

    if pszWhere is not None:
        poLyr.SetAttributeFilter(pszWhere)

    poFeat = poLyr.GetNextFeature()
    while poFeat is not None:
        poSrcGeom = poFeat.GetGeometryRef()
        if poSrcGeom is not None:
            eType = wkbFlatten(poSrcGeom.GetGeometryType())

            if poGeom is None:
                poGeom = ogr.Geometry( ogr.wkbMultiPolygon )

            if eType == ogr.wkbPolygon:
                poGeom.AddGeometry( poSrcGeom )
            elif eType == ogr.wkbMultiPolygon:
                for iGeom in range(poSrcGeom.GetGeometryCount()):
                    poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) )

            else:
                print("ERROR: Geometry not of polygon type." )
                if pszSQL is not None:
                    poDS.ReleaseResultSet( poLyr )
                poDS.Destroy()
                return None

        poFeat = poLyr.GetNextFeature()

    if pszSQL is not None:
        poDS.ReleaseResultSet( poLyr )
    poDS.Destroy()

    return poGeom

def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere):
    poGeom = None

    poDS = ogr.Open( pszDS, False )
    if poDS is None:
        return None

    if pszSQL is not None:
        poLyr = poDS.ExecuteSQL( pszSQL, None, None )
    elif pszLyr is not None:
        poLyr = poDS.GetLayerByName(pszLyr)
    else:
        poLyr = poDS.GetLayer(0)

    if poLyr is None:
        print("Failed to identify source layer from datasource.")
        poDS.Destroy()
        return None

    if pszWhere is not None:
        poLyr.SetAttributeFilter(pszWhere)

    poFeat = poLyr.GetNextFeature()
    while poFeat is not None:
        poSrcGeom = poFeat.GetGeometryRef()
        if poSrcGeom is not None:
            eType = wkbFlatten(poSrcGeom.GetGeometryType())

            if poGeom is None:
                poGeom = ogr.Geometry( ogr.wkbMultiPolygon )

            if eType == ogr.wkbPolygon:
                poGeom.AddGeometry( poSrcGeom )
            elif eType == ogr.wkbMultiPolygon:
                for iGeom in range(poSrcGeom.GetGeometryCount()):
                    poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) )

            else:
                print("ERROR: Geometry not of polygon type." )
                if pszSQL is not None:
                    poDS.ReleaseResultSet( poLyr )
                poDS.Destroy()
                return None

        poFeat = poLyr.GetNextFeature()

    if pszSQL is not None:
        poDS.ReleaseResultSet( poLyr )
    poDS.Destroy()

    return poGeom

def _create_layer(self, type, srid, attributes, title):
        """
        Creates an empty spatialite layer
        :param type: 'Point', 'LineString', 'Polygon', etc.
        :param srid: CRS ID of the layer
        :param attributes: list of (attribute_name, attribute_type, attribute_typename)
        :param title: title of the layer
        """
        driver = ogr.GetDriverByName('GPKG')
        ds = driver.CreateDataSource(self.output_local_file)
        layer = ds.CreateLayer("meta", geom_type = ogr.wkbNone)
        layer.CreateField(ogr.FieldDefn('key', ogr.OFTString))
        layer.CreateField(ogr.FieldDefn('value', ogr.OFTString))

        if srid:
            wkbType = { 'point': ogr.wkbPoint,
                        'multipoint': ogr.wkbMultiPoint,
                        'linestring': ogr.wkbLineString,
                        'multilinestring': ogr.wkbMultiLineString,
                        'polygon': ogr.wkbPolygon,
                        'multipolygon': ogr.wkbMultiPolygon
            }[type]
            srs = osr.SpatialReference()
            srs.ImportFromEPSGA(int(srid))
        else:
            wkbType = ogr.wkbNone
            srs = None
        layer = ds.CreateLayer("data", srs, wkbType, ['FID=id'])
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger64))
        layer.CreateField(ogr.FieldDefn('fid', ogr.OFTString))
        layer.CreateField(ogr.FieldDefn('_xml_', ogr.OFTString))

        att_type_map = {QVariant.String : ogr.OFTString,
                        QVariant.Int : ogr.OFTInteger,
                        QVariant.Double: ogr.OFTReal,
                        QVariant.DateTime: ogr.OFTDateTime}
        for aname, atype in attributes:
            layer.CreateField(ogr.FieldDefn(aname, att_type_map[atype]))

        # update fields
        layer.ResetReading()

        qgs_layer = QgsVectorLayer("{}|layername=data".format(self.output_local_file), title, "ogr")
        return qgs_layer

def split_shape_into_features(shape_name, destination_directory, column, name, sep):
    '''
    This method will take a input shape and iterate over its features, creating
    a new shape file with each one of them. It copies all the fields and the
    same spatial reference from the original file. The created files are saved
    in the destination directory using the number of the field given. 
    '''
    driver = ogr.GetDriverByName(str('ESRI Shapefile'))
    shape = driver.Open(shape_name, 0)
    layer = shape.GetLayer()
    layer_name = layer.GetName()
    spatial_reference = layer.GetSpatialRef()
    in_feature = layer.GetNextFeature()
    shape_files = []

    while in_feature:
        encoding = 'utf-8'
        in_feature_name = in_feature.GetField(column)
        in_name = create_filename_from_string(in_feature.GetField(name).decode(encoding))

        final_path = destination_directory + str(in_feature.GetField(0))
        create_directory_path(final_path)
        output_name = create_filename(final_path, '%s__%s.shp' % (in_feature_name, in_name))
        shape_files.append(output_name)

        if os.path.exists(output_name):
            driver.DeleteDataSource(output_name)
        datasource = driver.CreateDataSource(output_name)

        outLayer = datasource.CreateLayer(layer_name, spatial_reference, geom_type=ogr.wkbPolygon)

        inLayerDefn = layer.GetLayerDefn()
        for i in range(0, inLayerDefn.GetFieldCount()):
            fieldDefn = inLayerDefn.GetFieldDefn(i)
            #LOGGER.debug(fieldDefn.GetName())
            outLayer.CreateField(fieldDefn)

        outLayerDefn = outLayer.GetLayerDefn()
        geometry = in_feature.GetGeometryRef()

        out_feature = ogr.Feature(outLayerDefn)
        out_feature.SetGeometry(geometry)

        for i in range(0, outLayerDefn.GetFieldCount()):
            out_feature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), in_feature.GetField(i))

        outLayer.CreateFeature(out_feature)
        out_feature.Destroy()
        in_feature.Destroy()
        in_feature = layer.GetNextFeature()
    shape.Destroy()
    datasource.Destroy()
    return shape_files