Python ogr.GetDriverByName() 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 geotiff_dir():
    """

    Returns
    -------
    A pointer to a temporary folder that contains a 3-band geotiff
    of 3x3, with all values being 1.

    """
    tempDir = tempfile.TemporaryDirectory()
    fileformat = "GTiff"
    driver = gdal.GetDriverByName(fileformat)
    metadata = driver.GetMetadata()
    tempPath = os.path.join(tempDir.name)
    testDataset = driver.Create(os.path.join(tempDir.name, "tempTiff.tif"),
                                xsize=3, ysize=3, bands=3, eType=gdal.GDT_CFloat32)
    for i in range(3):
        testDataset.GetRasterBand(i+1).WriteArray(np.ones([3, 3]))
    testDataset = None
    yield tempPath
    tempDir.cleanup() 

def create_temp_tiff(self, name, content=np.ones([3, 3, 3]), geotransform=(10, 10, 0, 10, 0, -10)):
        """Creates a temporary geotiff in self.path
        """
        if len(content.shape) != 3:
            raise IndexError
        path = os.path.join(self.path, name)
        driver = gdal.GetDriverByName('GTiff')
        new_image = driver.Create(
            path,
            xsize=content.shape[1],
            ysize=content.shape[2],
            bands=content.shape[0],
            eType=gdal.GDT_Byte
        )
        new_image.SetGeoTransform(geotransform)
        for band in range(content.shape[0]):
            raster_band = new_image.GetRasterBand(band+1)
            raster_band.WriteArray(content[band, ...].T)
        new_image.SetProjection(self.srs.ExportToWkt())
        new_image.FlushCache()
        self.images.append(new_image)
        self.image_paths.append(path) 

def Get_Extent(input_lyr):
    """
    Obtain the input layer extent (xmin, ymin, xmax, ymax)
    """
    # Input
    filename, ext = os.path.splitext(input_lyr)
    if ext.lower() == '.shp':
        inp_driver = ogr.GetDriverByName('ESRI Shapefile')
        inp_source = inp_driver.Open(input_lyr)
        inp_lyr = inp_source.GetLayer()
        x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
        inp_lyr = None
        inp_source = None
    elif ext.lower() == '.tif':
        inp_lyr = gdal.Open(input_lyr)
        inp_transform = inp_lyr.GetGeoTransform()
        x_min = inp_transform[0]
        x_max = x_min + inp_transform[1] * inp_lyr.RasterXSize
        y_max = inp_transform[3]
        y_min = y_max + inp_transform[5] * inp_lyr.RasterYSize
        inp_lyr = None
    else:
        raise Exception('The input data type is not recognized')
    return (x_min, y_min, x_max, y_max) 

def Get_Extent(input_lyr):
    """
    Obtain the input layer extent (xmin, ymin, xmax, ymax)
    """
    # Input
    filename, ext = os.path.splitext(input_lyr)
    if ext.lower() == '.shp':
        inp_driver = ogr.GetDriverByName('ESRI Shapefile')
        inp_source = inp_driver.Open(input_lyr)
        inp_lyr = inp_source.GetLayer()
        x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
        inp_lyr = None
        inp_source = None
    elif ext.lower() == '.tif':
        inp_lyr = gdal.Open(input_lyr)
        inp_transform = inp_lyr.GetGeoTransform()
        x_min = inp_transform[0]
        x_max = x_min + inp_transform[1] * inp_lyr.RasterXSize
        y_max = inp_transform[3]
        y_min = y_max + inp_transform[5] * inp_lyr.RasterYSize
        inp_lyr = None
    else:
        raise Exception('The input data type is not recognized')
    return (x_min, y_min, x_max, y_max) 

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 List_Fields(input_lyr):
    """
    Lists the field names of input layer
    """
    # Input
    if isinstance(input_lyr, str):
        inp_driver = ogr.GetDriverByName('ESRI Shapefile')
        inp_source = inp_driver.Open(input_lyr, 0)
        inp_lyr = inp_source.GetLayer()
        inp_lyr_defn = inp_lyr.GetLayerDefn()
    elif isinstance(input_lyr, ogr.Layer):
        inp_lyr_defn = input_lyr.GetLayerDefn()

    # List
    names_ls = []

    # Loop
    for j in range(0, inp_lyr_defn.GetFieldCount()):
        field_defn = inp_lyr_defn.GetFieldDefn(j)
        names_ls.append(field_defn.GetName())

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

    # Return
    return names_ls 

def Array_to_Raster(input_array, output_tiff, ll_corner, cellsize,
                    srs_wkt):
    """
    Saves an array into a raster file
    """
    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    y_ncells, x_ncells = input_array.shape
    gdal_datatype = gdaltype_from_dtype(input_array.dtype)

    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal_datatype)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(-9999)

    out_top_left_x = ll_corner[0]
    out_top_left_y = ll_corner[1] + cellsize*y_ncells

    out_source.SetGeoTransform((out_top_left_x, cellsize, 0,
                                out_top_left_y, 0, -cellsize))
    out_source.SetProjection(str(srs_wkt))
    out_band.WriteArray(input_array)

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

    # Return
    return output_tiff 

def get_layer_names(geo_fpath, selected_layer_name=None):
    """
    Since layers in the same dataset can use
    different projections, we'll need to create a
    DbEntity for each layer. Use the osr library from
    gdal to loop through layers and grab their names.
    """

    _, geo_ext = os.path.splitext(geo_fpath)
    driver_name = OGR_DRIVER_MAP[geo_ext]

    driver = ogr.GetDriverByName(driver_name)
    dataset = driver.Open(geo_fpath, 0)

    layer_names = []
    for idx in xrange(dataset.GetLayerCount()):
        layer = dataset.GetLayerByIndex(idx)
        layer_name = layer.GetName()

        if selected_layer_name:
            if selected_layer_name == layer_name:
                layer_names.append(layer_name)
                break
        else:
            layer_names.append(layer_name)

    return layer_names 

def Extract_Band(input_tiff, output_tiff, band_number=1):
    """
    Extract and save a raster band into a new raster
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(band_number)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    NoData_value = inp_band.GetNoDataValue()

    x_ncells = inp_lyr.RasterXSize
    y_ncells = inp_lyr.RasterYSize

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform(inp_transform)
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(inp_array)

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

    # Return
    return output_tiff 

def Extract_Band(input_tiff, output_tiff, band_number=1):
    """
    Extract and save a raster band into a new raster
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(band_number)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    NoData_value = inp_band.GetNoDataValue()

    x_ncells = inp_lyr.RasterXSize
    y_ncells = inp_lyr.RasterYSize

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform(inp_transform)
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(inp_array)

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

    # Return
    return output_tiff 

def Array_to_Raster(input_array, output_tiff, ll_corner, cellsize,
                    srs_wkt):
    """
    Saves an array into a raster file
    """
    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    y_ncells, x_ncells = input_array.shape
    gdal_datatype = gdaltype_from_dtype(input_array.dtype)

    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal_datatype)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(-9999)

    out_top_left_x = ll_corner[0]
    out_top_left_y = ll_corner[1] + cellsize*y_ncells

    out_source.SetGeoTransform((out_top_left_x, cellsize, 0,
                                out_top_left_y, 0, -cellsize))
    out_source.SetProjection(str(srs_wkt))
    out_band.WriteArray(input_array)

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

    # Return
    return output_tiff 

def List_Fields(input_lyr):
    """
    Lists the field names of input layer
    """
    # Input
    if isinstance(input_lyr, str):
        inp_driver = ogr.GetDriverByName('ESRI Shapefile')
        inp_source = inp_driver.Open(input_lyr, 0)
        inp_lyr = inp_source.GetLayer()
        inp_lyr_defn = inp_lyr.GetLayerDefn()
    elif isinstance(input_lyr, ogr.Layer):
        inp_lyr_defn = input_lyr.GetLayerDefn()

    # List
    names_ls = []

    # Loop
    for j in range(0, inp_lyr_defn.GetFieldCount()):
        field_defn = inp_lyr_defn.GetFieldDefn(j)
        names_ls.append(field_defn.GetName())

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

    # Return
    return names_ls 

def Feature_to_Raster(input_shp, output_tiff, cellsize, field_name=False, NoData_value=-9999):
    """
    Converts a shapefile into a raster
    """

    # 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()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,1, gdal.GDT_Int16)
    print("+"*50)
    print(x_ncells, y_ncells,1, gdal.GDT_Int16)

    out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
    out_source.SetProjection(inp_srs.ExportToWkt())
    out_lyr = out_source.GetRasterBand(1)
    out_lyr.SetNoDataValue(NoData_value)

    # Rasterize
    # print(inp_lyr)
    if field_name:
        gdal.RasterizeLayer(out_source, [1], inp_lyr,options=["ATTRIBUTE={0}".format(field_name)])
    else:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])

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

    # Return
    return output_tiff 

#geo_silhouettes 

def saveToShape(array, srs, outShapeFile):
    # Parse a delimited text file of volcano data and create a shapefile
    # use a dictionary reader so we can access by field name
    # set up the shapefile driver
    import ogr
    outDriver = ogr.GetDriverByName('ESRI Shapefile')

    # create the data source
    if os.path.exists(outShapeFile):
        outDriver.DeleteDataSource(outShapeFile)
    # Remove output shapefile if it already exists

    # options = ['SPATIALITE=YES'])
    ds = outDriver.CreateDataSource(outShapeFile)

    # create the spatial reference, WGS84

    lyrout = ds.CreateLayer('randomSubset', srs, ogr.wkbPoint)
    fields = [array[1].GetFieldDefnRef(i).GetName()
              for i in range(array[1].GetFieldCount())]

    for i, f in enumerate(fields):
        field_name = ogr.FieldDefn(f, array[1].GetFieldDefnRef(i).GetType())
        field_name.SetWidth(array[1].GetFieldDefnRef(i).GetWidth())
        lyrout.CreateField(field_name)

    for k in array:
        lyrout.CreateFeature(k)

    # Save and close the data source
    ds = None 

def saveToShape(self, array, srs, outShapeFile):
        # Parse a delimited text file of volcano data and create a shapefile
        # use a dictionary reader so we can access by field name
        # set up the shapefile driver
        outDriver = ogr.GetDriverByName('ESRI Shapefile')

        # create the data source
        if os.path.exists(outShapeFile):
            outDriver.DeleteDataSource(outShapeFile)
        # Remove output shapefile if it already exists

        # options = ['SPATIALITE=YES'])
        ds = outDriver.CreateDataSource(outShapeFile)

        # create the spatial reference, WGS84

        lyrout = ds.CreateLayer('randomSubset', srs)
        fields = [
            array[1].GetFieldDefnRef(i).GetName() for i in range(
                array[1].GetFieldCount())]

        for f in fields:
            field_name = ogr.FieldDefn(f, ogr.OFTString)
            field_name.SetWidth(24)
            lyrout.CreateField(field_name)

        for k in array:
            lyrout.CreateFeature(k)

        # Save and close the data source
        ds = None 

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 Kriging_Interpolation_Points(input_shp, field_name, output_tiff, cellsize,
                                 bbox=None):
    """
    Interpolate point data using Ordinary Kriging

    Reference: https://cran.r-project.org/web/packages/automap/automap.pdf
    """
    # Spatial reference
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_srs = inp_lyr.GetSpatialRef()
    srs_wkt = inp_srs.ExportToWkt()
    inp_source = None
    # Temp folder
    temp_dir = tempfile.mkdtemp()
    temp_points_tiff = os.path.join(temp_dir, 'points_ras.tif')
    # Points to raster
    Feature_to_Raster(input_shp, temp_points_tiff,
                      cellsize, field_name, -9999)
    # Raster extent
    if bbox:
        xmin, ymin, xmax, ymax = bbox
        ll_corner = [xmin, ymin]
        x_ncells = int(math.ceil((xmax - xmin)/cellsize))
        y_ncells = int(math.ceil((ymax - ymin)/cellsize))
    else:
        temp_lyr = gdal.Open(temp_points_tiff)
        x_min, x_max, y_min, y_max = temp_lyr.GetExtent()
        ll_corner = [x_min, y_min]
        x_ncells = temp_lyr.RasterXSize
        y_ncells = temp_lyr.RasterYSize
        temp_lyr = None
    # Raster to array
    points_array = Raster_to_Array(temp_points_tiff, ll_corner,
                                   x_ncells, y_ncells, values_type='float32')
    # Run kriging
    x_vector = np.arange(xmin + cellsize/2, xmax + cellsize/2, cellsize)
    y_vector = np.arange(ymin + cellsize/2, ymax + cellsize/2, cellsize)
    out_array = Kriging_Interpolation_Array(points_array, x_vector, y_vector)
    # Save array as raster
    Array_to_Raster(out_array, output_tiff, ll_corner, cellsize, srs_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 Apply_Filter(input_tiff, output_tiff, number_of_passes):
    """
    Smooth a raster by replacing cell value by the average value of the
    surrounding cells
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    top_left_x = inp_transform[0]
    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    top_left_y = inp_transform[3]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    x_ncells = inp_lyr.RasterXSize
    y_ncells = inp_lyr.RasterYSize

    # Filter
    inp_array[inp_array == NoData_value] = pd.np.nan
    out_array = array_filter(inp_array, number_of_passes)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                top_left_y, rot_2, cellsize_y))
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(out_array)

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

    # Return
    return output_tiff 

def rasterRW(fn,raster_lyr,raster_lyr_w):
    gdal.UseExceptions()
    
    '''打开栅格数据'''
    try:
        src_ds=gdal.Open(os.path.join(fn,raster_lyr))
    except RuntimeError as e:
        print( 'Unable to open %s'% os.path.join(fn,raster_lyr))
        print(e)
        sys.exit(1)
    print("metadata:",src_ds.GetMetadata())   
    
    '''获取所有波段'''
    srcband=[]
    for band_num in range(1,5):
        try:
            srcband.append(src_ds.GetRasterBand(band_num))
        except RuntimeError as e:
            print('Band ( %i ) not found' % band_num)
            print(e)
            sys.exit(1)
    print(srcband)
    
    '''获取用于NDVI计算的红和近红波段数组,并计算ndvi'''
    red=srcband[0].ReadAsArray().astype(np.float)
    nir=srcband[3].ReadAsArray()
    red=np.ma.masked_where(nir+red==0,red)
    ndvi=(nir-red)/(nir+red)
    ndvi=ndvi.filled(-99)
    print(ndvi.shape,ndvi.std())
    
    '''初始化输出栅格'''
    driver=gdal.GetDriverByName('GTiff')
    out_raster=driver.Create(os.path.join(fn,raster_lyr_w),src_ds.RasterXSize,src_ds.RasterYSize,1,gdal.GDT_Float64)
    out_raster.SetProjection(src_ds.GetProjection()) #设置投影与参考栅格同
    out_raster.SetGeoTransform(src_ds.GetGeoTransform()) #配置地理转换与参考栅格同
    
    '''将数组传给栅格波段，为栅格值'''
    out_band=out_raster.GetRasterBand(1)
    out_band.WriteArray(ndvi)
    
    '''设置overview'''
    overviews = pb.compute_overview_levels(out_raster.GetRasterBand(1))
    out_raster.BuildOverviews('average', overviews)
    
    '''清理缓存与移除数据源'''
    out_band.FlushCache()
    out_band.ComputeStatistics(False)
    del src_ds,out_raster,out_band 

def Raster_to_Array(input_tiff, ll_corner, x_ncells, y_ncells,
                    values_type='float32'):
    """
    Loads a raster into a numpy array
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_data_type = inp_band.DataType

    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    ll_x = ll_corner[0]
    ll_y = ll_corner[1]

    top_left_x = ll_x
    top_left_y = ll_y - cellsize_y*y_ncells

    # Change start point
    temp_path = tempfile.mkdtemp()
    temp_driver = gdal.GetDriverByName('GTiff')
    temp_tiff = os.path.join(temp_path, os.path.basename(input_tiff))
    temp_source = temp_driver.Create(temp_tiff, x_ncells, y_ncells,
                                     1, inp_data_type)
    temp_source.GetRasterBand(1).SetNoDataValue(NoData_value)
    temp_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                 top_left_y, rot_2, cellsize_y))
    temp_source.SetProjection(inp_srs)

    # Snap
    gdal.ReprojectImage(inp_lyr, temp_source, inp_srs, inp_srs,
                        gdal.GRA_Bilinear)
    temp_source = None

    # Read array
    d_type = pd.np.dtype(values_type)
    out_lyr = gdal.Open(temp_tiff)
    array = out_lyr.ReadAsArray(0, 0, out_lyr.RasterXSize,
                                out_lyr.RasterYSize).astype(d_type)
    array[pd.np.isclose(array, NoData_value)] = pd.np.nan
    out_lyr = None

    return array 

def Feature_to_Raster(input_shp, output_tiff,
                      cellsize, field_name=False, NoData_value=-9999):
    """
    Converts a shapefile into a raster
    """

    # 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()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal.GDT_Int16)

    out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
    out_source.SetProjection(inp_srs.ExportToWkt())
    out_lyr = out_source.GetRasterBand(1)
    out_lyr.SetNoDataValue(NoData_value)

    # Rasterize
    if field_name:
        gdal.RasterizeLayer(out_source, [1], inp_lyr,
                            options=["ATTRIBUTE={0}".format(field_name)])
    else:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])

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

    # 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 Kriging_Interpolation_Points(input_shp, field_name, output_tiff, cellsize,
                                 bbox=None):
    """
    Interpolate point data using Ordinary Kriging

    Reference: https://cran.r-project.org/web/packages/automap/automap.pdf
    """
    # Spatial reference
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_srs = inp_lyr.GetSpatialRef()
    srs_wkt = inp_srs.ExportToWkt()
    inp_source = None
    # Temp folder
    temp_dir = tempfile.mkdtemp()
    temp_points_tiff = os.path.join(temp_dir, 'points_ras.tif')
    # Points to raster
    Feature_to_Raster(input_shp, temp_points_tiff,
                      cellsize, field_name, -9999)
    # Raster extent
    if bbox:
        xmin, ymin, xmax, ymax = bbox
        ll_corner = [xmin, ymin]
        x_ncells = int(math.ceil((xmax - xmin)/cellsize))
        y_ncells = int(math.ceil((ymax - ymin)/cellsize))
    else:
        temp_lyr = gdal.Open(temp_points_tiff)
        x_min, x_max, y_min, y_max = temp_lyr.GetExtent()
        ll_corner = [x_min, y_min]
        x_ncells = temp_lyr.RasterXSize
        y_ncells = temp_lyr.RasterYSize
        temp_lyr = None
    # Raster to array
    points_array = Raster_to_Array(temp_points_tiff, ll_corner,
                                   x_ncells, y_ncells, values_type='float32')
    # Run kriging
    x_vector = np.arange(xmin + cellsize/2, xmax + cellsize/2, cellsize)
    y_vector = np.arange(ymin + cellsize/2, ymax + cellsize/2, cellsize)
    out_array = Kriging_Interpolation_Array(points_array, x_vector, y_vector)
    # Save array as raster
    Array_to_Raster(out_array, output_tiff, ll_corner, cellsize, srs_wkt)
    # Return
    return output_tiff 

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 Apply_Filter(input_tiff, output_tiff, number_of_passes):
    """
    Smooth a raster by replacing cell value by the average value of the
    surrounding cells
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    top_left_x = inp_transform[0]
    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    top_left_y = inp_transform[3]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    x_ncells = inp_lyr.RasterXSize
    y_ncells = inp_lyr.RasterYSize

    # Filter
    inp_array[inp_array == NoData_value] = pd.np.nan
    out_array = array_filter(inp_array, number_of_passes)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                top_left_y, rot_2, cellsize_y))
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(out_array)

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

    # Return
    return output_tiff 

def Raster_to_Array(input_tiff, ll_corner, x_ncells, y_ncells,
                    values_type='float32'):
    """
    Loads a raster into a numpy array
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_data_type = inp_band.DataType

    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    ll_x = ll_corner[0]
    ll_y = ll_corner[1]

    top_left_x = ll_x
    top_left_y = ll_y - cellsize_y*y_ncells

    # Change start point
    temp_path = tempfile.mkdtemp()
    temp_driver = gdal.GetDriverByName('GTiff')
    temp_tiff = os.path.join(temp_path, os.path.basename(input_tiff))
    temp_source = temp_driver.Create(temp_tiff, x_ncells, y_ncells,
                                     1, inp_data_type)
    temp_source.GetRasterBand(1).SetNoDataValue(NoData_value)
    temp_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                 top_left_y, rot_2, cellsize_y))
    temp_source.SetProjection(inp_srs)

    # Snap
    gdal.ReprojectImage(inp_lyr, temp_source, inp_srs, inp_srs,
                        gdal.GRA_Bilinear)
    temp_source = None

    # Read array
    d_type = pd.np.dtype(values_type)
    out_lyr = gdal.Open(temp_tiff)
    array = out_lyr.ReadAsArray(0, 0, out_lyr.RasterXSize,
                                out_lyr.RasterYSize).astype(d_type)
    array[pd.np.isclose(array, NoData_value)] = pd.np.nan
    out_lyr = None

    return array 

def Feature_to_Raster(input_shp, output_tiff,
                      cellsize, field_name=False, NoData_value=-9999):
    """
    Converts a shapefile into a raster
    """

    # 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()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal.GDT_Int16)

    out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
    out_source.SetProjection(inp_srs.ExportToWkt())
    out_lyr = out_source.GetRasterBand(1)
    out_lyr.SetNoDataValue(NoData_value)

    # Rasterize
    if field_name:
        gdal.RasterizeLayer(out_source, [1], inp_lyr,
                            options=["ATTRIBUTE={0}".format(field_name)])
    else:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])

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

    # 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