Python ogr.Feature() Examples

def create_temp_shape(self, name, point_list):
        vector_file = os.path.join(self.temp_dir.name, name)
        shape_driver = ogr.GetDriverByName("ESRI Shapefile")  # Depreciated; replace at some point
        vector_data_source = shape_driver.CreateDataSource(vector_file)
        vector_layer = vector_data_source.CreateLayer("geometry", self.srs, geom_type=ogr.wkbPolygon)
        ring = ogr.Geometry(ogr.wkbLinearRing)
        for point in point_list:
            ring.AddPoint(point[0], point[1])
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        vector_feature_definition = vector_layer.GetLayerDefn()
        vector_feature = ogr.Feature(vector_feature_definition)
        vector_feature.SetGeometry(poly)
        vector_layer.CreateFeature(vector_feature)
        vector_layer.CreateField(ogr.FieldDefn("class", ogr.OFTInteger))
        feature = ogr.Feature(vector_layer.GetLayerDefn())
        feature.SetField("class", 3)

        vector_data_source.FlushCache()
        self.vectors.append(vector_data_source)  # Check this is the right thing to be saving here
        self.vector_paths.append(vector_file) 

def create_100x100_shp(self, name):
        """Cretes  a shapefile with a vector layer named "geometry" containing a 100mx100m square , top left corner
        being at wgs coords 10,10.
        This polygon has a field, 'class' with a value of 3. Left in for back-compatability"""
        # TODO Generalise this
        vector_file = os.path.join(self.temp_dir.name, name)
        shape_driver = ogr.GetDriverByName("ESRI Shapefile")  # Depreciated; replace at some point
        vector_data_source = shape_driver.CreateDataSource(vector_file)
        vector_layer = vector_data_source.CreateLayer("geometry", self.srs, geom_type=ogr.wkbPolygon)
        ring = ogr.Geometry(ogr.wkbLinearRing)
        ring.AddPoint(10.0, 10.0)
        ring.AddPoint(10.0, 110.0)
        ring.AddPoint(110.0, 110.0)
        ring.AddPoint(110.0, 10.0)
        ring.AddPoint(10.0, 10.0)
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        vector_feature_definition = vector_layer.GetLayerDefn()
        vector_feature = ogr.Feature(vector_feature_definition)
        vector_feature.SetGeometry(poly)
        vector_layer.CreateFeature(vector_feature)
        vector_layer.CreateField(ogr.FieldDefn("class", ogr.OFTInteger))
        feature = ogr.Feature(vector_layer.GetLayerDefn())
        feature.SetField("class", 3)

        vector_data_source.FlushCache()
        self.vectors.append(vector_data_source)  # Check this is the right thing to be saving here
        self.vector_paths.append(vector_file) 

def _map_column_values(self, feature, feature_defn, source_cols):
        """
        Retrieves values for specific columns from the specified feature.
        :param feature: Input feature.
        :type feature: ogr.Feature
        :param feature_defn: Feature definition for the layer.
        :type feature_defn: ogr.FeatureDefn
        :param source_cols: List of columns whose respective values will be
        retrieved.
        :type source_cols: list
        :return: Collection containing pairs of column names and corresponding
        values.
        :rtype: dict
        """
        col_values = {}

        if len(source_cols) == 0:
            return col_values

        for f in range(feature_defn.GetFieldCount()):
            field_defn = feature_defn.GetFieldDefn(f)
            field_name = field_defn.GetNameRef()

            match_idx = getIndex(source_cols, field_name)
            if match_idx != -1:
                field_value = feature.GetField(f)

                col_values[field_name] = field_value

        return col_values 

def save_point_list_to_shapefile(class_sample_point_dict, out_path, geotransform, projection_wkt, produce_csv=False):
    """Saves a list of points to a shapefile at out_path. Need the gt and projection of the raster.
    GT is needed to move each point to the centre of the pixel. Can also produce a .csv file for CoolEarth"""
    log = logging.getLogger(__name__)
    log.info("Saving point list to shapefile")
    log.debug("GT: {}\nProjection: {}".format(geotransform, projection_wkt))
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(out_path)
    srs = osr.SpatialReference()
    srs.ImportFromWkt(projection_wkt)
    layer = data_source.CreateLayer("validation_points", srs, ogr.wkbPoint)
    class_field = ogr.FieldDefn("class", ogr.OFTString)
    class_field.SetWidth(24)
    layer.CreateField(class_field)

    for map_class, point_list in class_sample_point_dict.items():
        for point in point_list:
            feature = ogr.Feature(layer.GetLayerDefn())
            coord = pyeo.coordinate_manipulation.pixel_to_point_coordinates(point, geotransform)
            offset = geotransform[1]/2   # Adds half a pixel offset so points end up in the center of pixels
            wkt = "POINT({} {})".format(coord[0]+offset, coord[1]-offset) # Never forget about negative y values in gts.
            new_point = ogr.CreateGeometryFromWkt(wkt)
            feature.SetGeometry(new_point)
            feature.SetField("class", map_class)
            layer.CreateFeature(feature)
            feature = None

    layer = None
    data_source = None

    if produce_csv:
        csv_out_path = out_path.rsplit('.')[0] + ".csv"
        with open(csv_out_path, "w") as csv_file:
            writer = csv.writer(csv_file)
            writer.writerow(["id", "yCoordinate", "xCoordinate"])

            # Join all points create single dimesional list of points (and revise the '*' operator)
            for id,  point in enumerate(itertools.chain(*class_sample_point_dict.values())):
                coord = pyeo.coordinate_manipulation.pixel_to_point_coordinates(point, geotransform)
                offset = geotransform[1] / 2  # Adds half a pixel offset so points end up in the center of pixels
                lat = coord[0] + offset
                lon = coord[1] - offset
                writer.writerow([id, lon, lat])
        log.info("CSV out at: {}".format(csv_out_path)) 

def Buffer(input_shp, output_shp, distance):
    """
    Creates a buffer of the input shapefile by a given distance
    """

    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_lyr_defn = inp_lyr.GetLayerDefn()
    inp_srs = inp_lyr.GetSpatialRef()

    # Output
    out_name = os.path.splitext(os.path.basename(output_shp))[0]
    out_driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(output_shp):
        out_driver.DeleteDataSource(output_shp)
    out_source = out_driver.CreateDataSource(output_shp)

    out_lyr = out_source.CreateLayer(out_name, inp_srs, ogr.wkbPolygon)
    out_lyr_defn = out_lyr.GetLayerDefn()

    # Add fields
    for i in range(inp_lyr_defn.GetFieldCount()):
        field_defn = inp_lyr_defn.GetFieldDefn(i)
        out_lyr.CreateField(field_defn)

    # Add features
    for i in range(inp_lyr.GetFeatureCount()):
        feature_inp = inp_lyr.GetNextFeature()
        geometry = feature_inp.geometry()
        feature_out = ogr.Feature(out_lyr_defn)

        for j in range(0, out_lyr_defn.GetFieldCount()):
            feature_out.SetField(out_lyr_defn.GetFieldDefn(j).GetNameRef(),
                                 feature_inp.GetField(j))

        feature_out.SetGeometry(geometry.Buffer(distance))
        out_lyr.CreateFeature(feature_out)
        feature_out = None

    # Save and/or close the data sources
    inp_source = None
    out_source = None

    # Return
    return output_shp 

def Interpolation_Default(input_shp, field_name, output_tiff,
                          method='nearest', cellsize=None):
    '''
    Interpolate point data into a raster

    Available methods: 'nearest', 'linear', 'cubic'
    '''
    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_srs = inp_lyr.GetSpatialRef()
    inp_wkt = inp_srs.ExportToWkt()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    ll_corner = [x_min, y_min]
    if not cellsize:
        cellsize = min(x_max - x_min, y_max - y_min)/25.0
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Feature points
    x = []
    y = []
    z = []
    for i in range(inp_lyr.GetFeatureCount()):
        feature_inp = inp_lyr.GetNextFeature()
        point_inp = feature_inp.geometry().GetPoint()

        x.append(point_inp[0])
        y.append(point_inp[1])
        z.append(feature_inp.GetField(field_name))

    x = pd.np.array(x)
    y = pd.np.array(y)
    z = pd.np.array(z)

    # Grid
    X, Y = pd.np.meshgrid(pd.np.linspace(x_min + cellsize/2.0,
                                         x_max - cellsize/2.0,
                                         x_ncells),
                          pd.np.linspace(y_min + cellsize/2.0,
                                         y_max - cellsize/2.0,
                                         y_ncells))

    # Interpolate
    out_array = griddata((x, y), z, (X, Y), method=method)
    out_array = pd.np.flipud(out_array)

    # Save raster
    Array_to_Raster(out_array, output_tiff, ll_corner, cellsize, inp_wkt)

    # Return
    return output_tiff 

def Buffer(input_shp, output_shp, distance):
    """
    Creates a buffer of the input shapefile by a given distance
    """

    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_lyr_defn = inp_lyr.GetLayerDefn()
    inp_srs = inp_lyr.GetSpatialRef()

    # Output
    out_name = os.path.splitext(os.path.basename(output_shp))[0]
    out_driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(output_shp):
        out_driver.DeleteDataSource(output_shp)
    out_source = out_driver.CreateDataSource(output_shp)

    out_lyr = out_source.CreateLayer(out_name, inp_srs, ogr.wkbPolygon)
    out_lyr_defn = out_lyr.GetLayerDefn()

    # Add fields
    for i in range(inp_lyr_defn.GetFieldCount()):
        field_defn = inp_lyr_defn.GetFieldDefn(i)
        out_lyr.CreateField(field_defn)

    # Add features
    for i in range(inp_lyr.GetFeatureCount()):
        feature_inp = inp_lyr.GetNextFeature()
        geometry = feature_inp.geometry()
        feature_out = ogr.Feature(out_lyr_defn)

        for j in range(0, out_lyr_defn.GetFieldCount()):
            feature_out.SetField(out_lyr_defn.GetFieldDefn(j).GetNameRef(),
                                 feature_inp.GetField(j))

        feature_out.SetGeometry(geometry.Buffer(distance))
        out_lyr.CreateFeature(feature_out)
        feature_out = None

    # Save and/or close the data sources
    inp_source = None
    out_source = None

    # Return
    return output_shp 

def Interpolation_Default(input_shp, field_name, output_tiff,
                          method='nearest', cellsize=None):
    '''
    Interpolate point data into a raster

    Available methods: 'nearest', 'linear', 'cubic'
    '''
    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_srs = inp_lyr.GetSpatialRef()
    inp_wkt = inp_srs.ExportToWkt()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    ll_corner = [x_min, y_min]
    if not cellsize:
        cellsize = min(x_max - x_min, y_max - y_min)/25.0
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Feature points
    x = []
    y = []
    z = []
    for i in range(inp_lyr.GetFeatureCount()):
        feature_inp = inp_lyr.GetNextFeature()
        point_inp = feature_inp.geometry().GetPoint()

        x.append(point_inp[0])
        y.append(point_inp[1])
        z.append(feature_inp.GetField(field_name))

    x = pd.np.array(x)
    y = pd.np.array(y)
    z = pd.np.array(z)

    # Grid
    X, Y = pd.np.meshgrid(pd.np.linspace(x_min + cellsize/2.0,
                                         x_max - cellsize/2.0,
                                         x_ncells),
                          pd.np.linspace(y_min + cellsize/2.0,
                                         y_max - cellsize/2.0,
                                         y_ncells))

    # Interpolate
    out_array = griddata((x, y), z, (X, Y), method=method)
    out_array = pd.np.flipud(out_array)

    # Save raster
    Array_to_Raster(out_array, output_tiff, ll_corner, cellsize, inp_wkt)

    # Return
    return output_tiff 

def pointWriting(fn,pt_lyrName_w,ref_lyr=False):
    ds=ogr.Open(fn,1)
    
    '''参考层，用于空间坐标投影，字段属性等参照'''
    ref_lyr=ds.GetLayer(ref_lyr)
    ref_sr=ref_lyr.GetSpatialRef()
    print(ref_sr)
    ref_schema=ref_lyr.schema #查看属性表字段名和类型
    for field in ref_schema:
        print(field.name,field.GetTypeName())   
 
    '''建立新的datasource数据源'''
    sf_driver=ogr.GetDriverByName('ESRI Shapefile')
    sfDS=os.path.join(fn,r'sf')
    if os.path.exists(sfDS):
        sf_driver.DeleteDataSource(sfDS)
    pt_ds=sf_driver.CreateDataSource(sfDS)
    if pt_ds is None:
        sys.exit('Could not open{0}'.format(sfDS))
        
    '''建立新layer层'''    
    if pt_ds.GetLayer(pt_lyrName_w):
        pt_ds.DeleteLayer(pt_lyrName_w)    
    pt_lyr=pt_ds.CreateLayer(pt_lyrName_w,ref_sr,ogr.wkbPoint)
    
    '''配置字段，名称以及类型和相关参数'''
    pt_lyr.CreateFields(ref_schema)
    LatFd=ogr.FieldDefn("origiLat",ogr.OFTReal)
    LatFd.SetWidth(8)
    LatFd.SetPrecision(3)
    pt_lyr.CreateField(LatFd)
    LatFd.SetName("Lat")
    pt_lyr.CreateField(LatFd)
     
    '''建立feature空特征和设置geometry几何类型'''
    print(pt_lyr.GetLayerDefn())
    pt_feat=ogr.Feature(pt_lyr.GetLayerDefn())    
    
    for feat in ref_lyr:  #循环feature
        '''设置几何体'''
        pt_ref=feat.geometry().Clone()
        wkt="POINT(%f %f)" %  (float(pt_ref.GetX()+0.01) , float(pt_ref.GetY()+0.01))
        newPt=ogr.CreateGeometryFromWkt(wkt) #使用wkt的方法建立点
        pt_feat.SetGeometry(newPt)
        '''设置字段值'''
        for i_field in range(feat.GetFieldCount()):
            pt_feat.SetField(i_field,feat.GetField(i_field))
        pt_feat.SetField("origiLat",pt_ref.GetX())
        pt_feat.SetField("Lat",pt_ref.GetX()+0.01)
        '''根据设置的几何体和字段值，建立feature。循环建立多个feature特征'''
        pt_lyr.CreateFeature(pt_feat)    
    del ds 

def lineWriting(fn,ln_lyrName_w,ref_lyr=False):
    ds=ogr.Open(fn,1)
    
    '''参考层，用于空间坐标投影，字段属性等参照'''
    ref_lyr=ds.GetLayer(ref_lyr)
    ref_sr=ref_lyr.GetSpatialRef()
    print(ref_sr)
    ref_schema=ref_lyr.schema #查看属性表字段名和类型
    for field in ref_schema:
        print(field.name,field.GetTypeName())      

    '''建立新layer层'''    
    if ds.GetLayer(ln_lyrName_w):
        ds.DeleteLayer(ln_lyrName_w)    
    ln_lyr=ds.CreateLayer(ln_lyrName_w,ref_sr,ogr.wkbMultiLineString)    

    '''配置字段，名称以及类型和相关参数'''
    Fd=ogr.FieldDefn("length",ogr.OFTReal)
    Fd.SetWidth(8)
    Fd.SetPrecision(3)
    ln_lyr.CreateField(Fd)
    Fd=ogr.FieldDefn("name",ogr.OFTString)
    ln_lyr.CreateField(Fd)    

    '''建立feature空特征和设置geometry几何类型'''
    print(ln_lyr.GetLayerDefn())
    ln_feat=ogr.Feature(ln_lyr.GetLayerDefn())

    for feat in ref_lyr:  #循环feature
        '''设置几何体'''
        ln_ref=feat.geometry().Clone()
        temp=""
        for j in range(ln_ref.GetPointCount()):
            if j==ln_ref.GetPointCount()-1:
                temp+="%f %f"%(float(ln_ref.GetX(j)+0.01) , float(ln_ref.GetY(j)+0.01))
            else:
                temp+="%f %f,"%(float(ln_ref.GetX(j)+0.01) , float(ln_ref.GetY(j)+0.01))
        wkt="LINESTRING(%s)" % temp  #使用wkt的方法建立直线
#        print(wkt)
        newLn=ogr.CreateGeometryFromWkt(wkt)        
        ln_feat.SetGeometry(newLn)
        '''设置字段值'''
        ln_feat.SetField("name",feat.GetField("name"))
        ln_feat.SetField("length",newLn.Length())
        '''根据设置的几何体和字段值，建立feature。循环建立多个feature特征'''
        ln_lyr.CreateFeature(ln_feat)    
    del ds 

def polygonWriting(fn,pg_lyrName_w,ref_lyr=False):
    ds=ogr.Open(fn,1)
    
    '''参考层，用于空间坐标投影，字段属性等参照'''
    ref_lyr=ds.GetLayer(ref_lyr)
    ref_sr=ref_lyr.GetSpatialRef()
#    print(ref_sr)
    ref_schema=ref_lyr.schema #查看属性表字段名和类型
    for field in ref_schema:
        print(field.name,field.GetTypeName())      

    '''建立新layer层'''    
    if ds.GetLayer(pg_lyrName_w):
        ds.DeleteLayer(pg_lyrName_w)    
    pg_lyr=ds.CreateLayer(pg_lyrName_w,ref_sr,ogr.wkbMultiPolygon)    

    '''配置字段，名称以及类型和相关参数'''
    Fd=ogr.FieldDefn("area",ogr.OFTReal)
    Fd.SetWidth(8)
    Fd.SetPrecision(8)
    pg_lyr.CreateField(Fd)
    Fd=ogr.FieldDefn("name",ogr.OFTString)
    pg_lyr.CreateField(Fd)    

    '''建立feature空特征和设置geometry几何类型'''
#    print(pg_lyr.GetLayerDefn())
    pg_feat=ogr.Feature(pg_lyr.GetLayerDefn())
    
    for feat in ref_lyr:  #循环feature
        '''设置几何体'''
        pg_ref=feat.geometry().Clone()
        tempSub=""
        for j in range(pg_ref.GetGeometryCount()):
            ring=pg_ref.GetGeometryRef(j)
            vertexes=ring.GetPoints()
#            print(len(vertexes))
            temp=""
            for i in range(len(vertexes)):                
                if i==len(vertexes)-1:
                    temp+="%f %f"%(float(vertexes[i][0]+0.01) , float(vertexes[i][1]+0.01))
                else:
                    temp+="%f %f,"%(float(vertexes[i][0]+0.01) , float(vertexes[i][1]+0.01))
            if j==pg_ref.GetGeometryCount()-1:
                tempSub+="(%s)"%temp   
            else:
                tempSub+="(%s),"%temp
#        print(tempSub)    
        wkt="POLYGON(%s)" % tempSub  #使用wkt的方法建立直线
#        print(wkt)
        newPg=ogr.CreateGeometryFromWkt(wkt)        
        pg_feat.SetGeometry(newPg)
        
        '''设置字段值'''
        pg_feat.SetField("name",feat.GetField("NAME"))
        pg_feat.SetField("area",newPg.GetArea())
        '''根据设置的几何体和字段值，建立feature。循环建立多个feature特征'''
        pg_lyr.CreateFeature(pg_feat)    
    del ds