Python osgeo.gdal 模块，GetDriverByName() 实例源码

def createGeoJSONFromRaster(geoJsonFileName, array2d, geom, proj,
                            layerName="BuildingID",
                            fieldName="BuildingID"):

    memdrv = gdal.GetDriverByName('MEM')
    src_ds = memdrv.Create('', array2d.shape[1], array2d.shape[0], 1)
    src_ds.SetGeoTransform(geom)
    src_ds.SetProjection(proj)
    band = src_ds.GetRasterBand(1)
    band.WriteArray(array2d)

    dst_layername = "BuildingID"
    drv = ogr.GetDriverByName("geojson")
    dst_ds = drv.CreateDataSource(geoJsonFileName)
    dst_layer = dst_ds.CreateLayer(layerName, srs=None)

    fd = ogr.FieldDefn(fieldName, ogr.OFTInteger)
    dst_layer.CreateField(fd)
    dst_field = 1

    gdal.Polygonize(band, None, dst_layer, dst_field, [], callback=None)

    return

def reprojectRaster(src,match_ds,dst_filename):
    src_proj = src.GetProjection()
    src_geotrans = src.GetGeoTransform()

    match_proj = match_ds.GetProjection()
    match_geotrans = match_ds.GetGeoTransform()
    wide = match_ds.RasterXSize
    high = match_ds.RasterYSize

    dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Int16)
    dst.SetGeoTransform( match_geotrans )
    dst.SetProjection( match_proj)

    gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)

    del dst # Flush
    return(gdal.Open(dst_filename,gdalconst.GA_ReadOnly))


#finds the intersecting extent of a series of scenes (left,right,bottom,top are arrays containing the respective lat/lon of the left,right,bottom,top of each image)

def createRasterFromGeoJson(srcGeoJson, srcRasterFileName, outRasterFileName):
    NoData_value = 0
    source_ds = ogr.Open(srcGeoJson)
    source_layer = source_ds.GetLayer()

    srcRaster = gdal.Open(srcRasterFileName)


    # Create the destination data source
    target_ds = gdal.GetDriverByName('GTiff').Create(outRasterFileName, srcRaster.RasterXSize, srcRaster.RasterYSize, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform(srcRaster.GetGeoTransform())
    target_ds.SetProjection(srcRaster.GetProjection())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(NoData_value)

    # Rasterize
    gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])
    band.FlushCache()

def reproject_tiff_from_model(old_name, new_name, model, unit):
    """
    Reprojects an tiff on a tiff model. Can be used to warp tiff.
    Keyword arguments:
    old_name -- the name of the old tiff file
    new_name -- the name of the output tiff file
    model -- the gdal dataset which will be used to warp the tiff
    unit -- the gdal unit in which the operation will be performed
    """
    mem_drv = gdal.GetDriverByName("MEM")

    old = gdal.Open(old_name)

    new = mem_drv.Create(new_name, model.RasterXSize, model.RasterYSize, 1,
                         unit)

    new.SetGeoTransform(model.GetGeoTransform())
    new.SetProjection(model.GetProjection())

    res = gdal.ReprojectImage(old, new, old.GetProjection(),
                              model.GetProjection(), gdal.GRA_NearestNeighbour)

    assert res == 0, 'Error in ReprojectImage'
    arr = new.ReadAsArray()
    new = None
    return arr

def write_int16_to_tiff(name, data, sr, geot, nodata_val=None):
    assert data.dtype == np.int16
    gtiff_drv = gdal.GetDriverByName("GTiff")
    tiff_file = gtiff_drv.Create(name, data.shape[1], data.shape[0], 1,
                                 gdal.GDT_Int16,
                                 options=['COMPRESS=DEFLATE', 'ZLEVEL=1'])
    tiff_file.SetGeoTransform(geot)
    tiff_file.SetProjection(sr)

    band = tiff_file.GetRasterBand(1)
    if nodata_val is not None:
        band.SetNoDataValue(nodata_val)
    band.WriteArray(data)
    band.FlushCache()
    del band
    del tiff_file

def create_mask_from_vector(vector_data_path, cols, rows, geo_transform, projection, target_value=1,
                            output_fname='', dataset_format='MEM'):
    """
    Rasterize the given vector (wrapper for gdal.RasterizeLayer). Return a gdal.Dataset.
    :param vector_data_path: Path to a shapefile
    :param cols: Number of columns of the result
    :param rows: Number of rows of the result
    :param geo_transform: Returned value of gdal.Dataset.GetGeoTransform (coefficients for
                          transforming between pixel/line (P,L) raster space, and projection
                          coordinates (Xp,Yp) space.
    :param projection: Projection definition string (Returned by gdal.Dataset.GetProjectionRef)
    :param target_value: Pixel value for the pixels. Must be a valid gdal.GDT_UInt16 value.
    :param output_fname: If the dataset_format is GeoTIFF, this is the output file name
    :param dataset_format: The gdal.Dataset driver name. [default: MEM]
    """
    data_source = gdal.OpenEx(vector_data_path, gdal.OF_VECTOR)
    if data_source is None:
        report_and_exit("File read failed: %s", vector_data_path)
    layer = data_source.GetLayer(0)
    driver = gdal.GetDriverByName(dataset_format)
    target_ds = driver.Create(output_fname, cols, rows, 1, gdal.GDT_UInt16)
    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(projection)
    gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[target_value])
    return target_ds

def write_stats(self):
        #if not hasattr(self, 'stack_count'):
        stat_list = ['stack_count', 'stack_mean', 'stack_std', 'stack_min', 'stack_max']
        if self.med:
            stat_list.append('stack_med')
            stat_list.append('stack_nmad')
        if any([not hasattr(self, i) for i in stat_list]):
            self.compute_stats()
        print("Writing out stats")
        #Create dummy ds - might want to use vrt here instead
        driver = gdal.GetDriverByName("MEM")
        ds = driver.Create('', self.ma_stack.shape[2], self.ma_stack.shape[1], 1, gdal.GDT_Float32)
        ds.SetGeoTransform(self.gt)
        ds.SetProjection(self.proj)
        #Write out with malib, should preserve ma type
        out_prefix = os.path.splitext(self.stack_fn)[0]
        iolib.writeGTiff(self.stack_count, out_prefix+'_count.tif', ds)
        iolib.writeGTiff(self.stack_mean, out_prefix+'_mean.tif', ds)
        iolib.writeGTiff(self.stack_std, out_prefix+'_std.tif', ds)
        iolib.writeGTiff(self.stack_min, out_prefix+'_min.tif', ds)
        iolib.writeGTiff(self.stack_max, out_prefix+'_max.tif', ds)
        if self.med:
            iolib.writeGTiff(self.stack_med, out_prefix+'_med.tif', ds)
            iolib.writeGTiff(self.stack_nmad, out_prefix+'_nmad.tif', ds)

def write_trend(self):
        #stat_list = ['stack_trend', 'stack_intercept', 'stack_detrended_std', 'stack_rsquared']
        stat_list = ['stack_trend', 'stack_intercept', 'stack_detrended_std']
        if any([not hasattr(self, i) for i in stat_list]):
            self.compute_trend()
        print("Writing out trend")
        #Create dummy ds - might want to use vrt here instead
        driver = gdal.GetDriverByName("MEM")
        ds = driver.Create('', self.ma_stack.shape[2], self.ma_stack.shape[1], 1, gdal.GDT_Float32)
        ds.SetGeoTransform(self.gt)
        ds.SetProjection(self.proj)
        #Write out with malib, should preserve ma type
        out_prefix = os.path.splitext(self.stack_fn)[0]
        iolib.writeGTiff(self.stack_trend, out_prefix+'_trend.tif', ds)
        iolib.writeGTiff(self.stack_intercept, out_prefix+'_intercept.tif', ds)
        iolib.writeGTiff(self.stack_detrended_std, out_prefix+'_detrended_std.tif', ds)
        #iolib.writeGTiff(self.stack_rsquared, out_prefix+'_rsquared.tif', ds)

def mem_ds(res, extent, srs=None, dtype=gdal.GDT_Float32):
    """Create a new GDAL Dataset in memory

    Useful for various applications that require a Dataset
    """
    #These round down to int
    #dst_ns = int((extent[2] - extent[0])/res)
    #dst_nl = int((extent[3] - extent[1])/res)
    #This should pad by 1 pixel, but not if extent and res were calculated together to give whole int
    dst_ns = int((extent[2] - extent[0])/res + 0.99)
    dst_nl = int((extent[3] - extent[1])/res + 0.99)
    m_ds = gdal.GetDriverByName('MEM').Create('', dst_ns, dst_nl, 1, dtype)
    m_gt = [extent[0], res, 0, extent[3], 0, -res]
    m_ds.SetGeoTransform(m_gt)
    if srs is not None:
        m_ds.SetProjection(srs.ExportToWkt())
    return m_ds

#Modify proj/gt of dst_fn in place

def array2Raster(result_arr,affine_par,dstFilename):
    # save arr to geotiff
    driver = gdal.GetDriverByName('GTiff')

    [A, D, B, E, C, F] = affine_par

    if result_arr.ndim == 2:
        [height,width] = result_arr.shape
        dataset = driver.Create(dstFilename, width, height, numBand, gdal.GDT_Float32)
        dataset.SetGeoTransform((C, A, B, F, D, E))
        dataset.GetRasterBand(1).WriteArray(array[:, :])
        #dataset.GetRasterBand(1).SetNoDataValue(0)

    elif result_arr.ndim == 3:
        [height,width,numBand] = result_arr.shape
        dataset = driver.Create(dstFilename, width, height, numBand, gdal.GDT_Float32)
        dataset.SetGeoTransform((C, A, B, F, D, E))
        for i in range(numBand):
            dataset.GetRasterBand(i + 1).WriteArray(result_arr[:, :, i])
            #dataset.GetRasterBand(i + 1).SetNoDataValue(-999.0)

    dataset.FlushCache()        # Write to disk

def rasterize(self,inRaster,inShape,inField):
        filename = tempfile.mktemp('.tif')        
        data = gdal.Open(inRaster,gdal.GA_ReadOnly)
        shp = ogr.Open(inShape)

        lyr = shp.GetLayer()

        driver = gdal.GetDriverByName('GTiff')
        dst_ds = driver.Create(filename,data.RasterXSize,data.RasterYSize, 1,gdal.GDT_Byte)
        dst_ds.SetGeoTransform(data.GetGeoTransform())
        dst_ds.SetProjection(data.GetProjection())
        OPTIONS = 'ATTRIBUTE='+inField
        gdal.RasterizeLayer(dst_ds, [1], lyr, None,options=[OPTIONS])
        data,dst_ds,shp,lyr=None,None,None,None


        return filename

def createImageDS(filename, x_min, y_min, x_max, y_max, pixel_size,  srs=None):
    # Create the destination data source
    x_res = int((x_max - x_min) / pixel_size) # resolution
    y_res = int((y_max - y_min) / pixel_size) # resolution
    ds = gdal.GetDriverByName('GTiff').Create(filename, x_res,
            y_res, 1, gdal.GDT_Byte)
    ds.SetGeoTransform((
            x_min, pixel_size, 0,
            y_max, 0, -pixel_size,
        ))
    if srs:
        # Make the target raster have the same projection as the source
        ds.SetProjection(srs.ExportToWkt())
    else:
        # Source has no projection (needs GDAL >= 1.7.0 to work)
        ds.SetProjection('LOCAL_CS["arbitrary"]')

    # Set nodata
    band = ds.GetRasterBand(1)
    band.SetNoDataValue(0)

    return ds

def rasterize_polygons(polygon_shp, template_raster, out_raster, field):
    """ generate a categorical raster based on polygons

    :rtype: None
    :param polygon_shp: input polygon shapefile
    :param template_raster: raster template for cellsize and extent
    :param out_raster: output raster file
    """

    # Source: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

    gdal.UseExceptions()
    # Get template raster specs
    src_ds = gdal.Open(template_raster)
    if src_ds is None:
        print 'Unable to open %s' % template_raster
        sys.exit(1)
    try:
        srcband = src_ds.GetRasterBand(1)
    except RuntimeError, e:
        print 'No band %i found' % 0
        print e
        sys.exit(1)

    # Open the data source and read in the extent
    source_ds = ogr.Open(polygon_shp)
    source_layer = source_ds.GetLayer()

    target_ds = gdal.GetDriverByName('GTiff').Create(out_raster, src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Float32)
    target_ds.SetGeoTransform(src_ds.GetGeoTransform())
    target_ds.SetProjection(src_ds.GetProjection())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(srcband.GetNoDataValue())

    # Rasterize
    gdal.RasterizeLayer(target_ds, [1], source_layer, options=["ATTRIBUTE={}".format(field)])

def array2raster(newRasterFilename,rasterOrigin,pixelWidth,pixelHeight,array,dataType):
    cols=array.shape[1]
    rows=array.shape[0]
    originX=rasterOrigin[0]
    originY=rasterOrigin[1]
    driver=gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterFilename, cols, rows, 1, dataType)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(32645)# this is the EPSG code for Nepal, should be changed for other locations
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

#takes a series of rasters and clips them to minExtent, then returns them as numpy arrays

def array2raster(newRasterFilename,rasterOrigin,pixelWidth,pixelHeight,array,dataType):
    array=array.astype(float)
    cols=array.shape[1]
    rows=array.shape[0]
    originX=rasterOrigin[0]
    originY=rasterOrigin[1]
    driver=gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterFilename, cols, rows, 1, dataType)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)#EPSG code for Nepal only
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()    

#gets the lat/lon extent of a raster

def polygonize(self,rasterTemp,outShp):

            sourceRaster = gdal.Open(rasterTemp)
            band = sourceRaster.GetRasterBand(1)
            driver = ogr.GetDriverByName("ESRI Shapefile")
            # If shapefile already exist, delete it
            if os.path.exists(outShp):
                driver.DeleteDataSource(outShp)

            outDatasource = driver.CreateDataSource(outShp)            
            # get proj from raster            
            srs = osr.SpatialReference()
            srs.ImportFromWkt( sourceRaster.GetProjectionRef() )
            # create layer with proj
            outLayer = outDatasource.CreateLayer(outShp,srs)
            # Add class column (1,2...) to shapefile

            newField = ogr.FieldDefn('Class', ogr.OFTInteger)
            outLayer.CreateField(newField)

            gdal.Polygonize(band, None,outLayer, 0,[],callback=None)  

            outDatasource.Destroy()
            sourceRaster=None
            band=None


            ioShpFile = ogr.Open(outShp, update = 1)


            lyr = ioShpFile.GetLayerByIndex(0)

            lyr.ResetReading()    

            for i in lyr:
                lyr.SetFeature(i)
            # if area is less than inMinSize or if it isn't forest, remove polygon 
                if i.GetField('Class')!=1:
                    lyr.DeleteFeature(i.GetFID())        
            ioShpFile.Destroy()

            return outShp

def average(input_dir, out_dir, size):
    input_dir = abspath(input_dir)
    log.info('Reading tifs from {}'.format(input_dir))
    tifs = glob.glob(join(input_dir, '*.tif'))

    for tif in tifs:
        data_source = gdal.Open(tif, gdal.GA_ReadOnly)
        band = data_source.GetRasterBand(1)
        # data_type = gdal.GetDataTypeName(band.DataType)
        data = band.ReadAsArray()
        no_data_val = band.GetNoDataValue()
        averaged_data = filter_data(data, size, no_data_val)
        log.info('Calculated average for {}'.format(basename(tif)))

        output_file = join(out_dir, 'average_' + basename(tif))
        out_ds = gdal.GetDriverByName('GTiff').Create(
            output_file, data_source.RasterXSize, data_source.RasterYSize,
            1, band.DataType)
        out_band = out_ds.GetRasterBand(1)
        out_band.WriteArray(averaged_data)
        out_ds.SetGeoTransform(data_source.GetGeoTransform())
        out_ds.SetProjection(data_source.GetProjection())
        out_band.FlushCache()  # Write data to disc
        out_ds = None  # close out_ds
        data_source = None  # close dataset

        log.info('Finished converting {}'.format(basename(tif)))

def write_file(output_name,output_array,GeoT,xsize,ysize,proJ,driverName='GTiff'):
    print "creating", output_name
    dr=gdal.GetDriverByName(driverName)
    dr.Register()
    do=dr.Create(output_name,xsize,ysize,1,gdal.GDT_UInt16,options = [ 'COMPRESS=LZW' ])

    do.SetGeoTransform(GeoT)
    do.SetProjection(proJ)
    do.GetRasterBand(1).WriteArray(output_array)
    do.GetRasterBand(1).SetNoDataValue(65535)
    do=None


# Function to wirte multi-dimesion array to tiff file

def write_file(output_name, output_array, GeoT, xsize, ysize, proJ, driverName='GTiff'):
    print "creating", output_name
    dr = gdal.GetDriverByName(driverName)
    dr.Register()
    do = dr.Create(output_name, xsize, ysize, 1, gdal.GDT_UInt16, options=['COMPRESS=LZW'])

    do.SetGeoTransform(GeoT)
    do.SetProjection(proJ)
    do.GetRasterBand(1).WriteArray(output_array)
    do.GetRasterBand(1).SetNoDataValue(32767)
    do = None

def write_file(output_name, output_array, GeoT, xsize, ysize, proJ, driverName='GTiff'):
    print "creating", output_name
    dr = gdal.GetDriverByName(driverName)
    dr.Register()
    do = dr.Create(output_name, xsize, ysize, 1, gdal.GDT_UInt16, options=['COMPRESS=LZW'])

    do.SetGeoTransform(GeoT)
    do.SetProjection(proJ)
    do.GetRasterBand(1).WriteArray(output_array)
    do.GetRasterBand(1).SetNoDataValue(65535)
    do = None

def write_file(output_name,output_array,GeoT,xsize,ysize,proJ,driverName='GTiff'):
    print "creating", output_name
    dr=gdal.GetDriverByName(driverName)
    dr.Register()
    do=dr.Create(output_name,xsize,ysize,1,gdal.GDT_Float32)
    do.SetGeoTransform(GeoT)
    do.SetProjection(proJ)
    do.GetRasterBand(1).WriteArray(output_array)
    do.GetRasterBand(1).SetNoDataValue(32767)
    do=None

def write_file(output_name,output_array,GeoT,xsize,ysize,proJ,driverName='GTiff'):
    print "creating", output_name
    dr=gdal.GetDriverByName(driverName)
    dr.Register()
    do=dr.Create(output_name,xsize,ysize,1,gdal.GDT_Float32)
    do.SetGeoTransform(GeoT)
    do.SetProjection(proJ)
    do.GetRasterBand(1).WriteArray(output_array)
    do.GetRasterBand(1).SetNoDataValue(32767)
    do=None

# Function to calculate the moving average

def write_file(output_name,output_array,GeoT,xsize,ysize,proJ,driverName='GTiff'):
    print "creating", output_name
    dr=gdal.GetDriverByName(driverName)
    dr.Register()
    do=dr.Create(output_name,xsize,ysize,1,gdal.GDT_Float32)
    do.SetGeoTransform(GeoT)
    do.SetProjection(proJ)
    do.GetRasterBand(1).WriteArray(output_array)
    do.GetRasterBand(1).SetNoDataValue(32767)
    do=None

def create_gpkg(
    gpkg_name, proj_string, size=(1, 1), geotransform=[0, 1, 0, 0, 0, -1],
    creation_options=None
):
    if os.path.exists("%s.gpkg" % gpkg_name):
        sys.stderr.write(
            "ERROR: SQLite GeoPackage '%s.gpkg' already exists.\n" % gpkg_name
        )
        sys.exit(1)

    gdal.AllRegister()
    drv = gdal.GetDriverByName("GPKG")
    try:
        gpkg = drv.Create(
            "%s.gpkg" % gpkg_name, size[0], size[1], 1, gdal.GDT_Byte,
            creation_options
        )

        proj = osr.SpatialReference()
        res = proj.SetWellKnownGeogCS(proj_string)
        if res != 0:
            if proj_string[0:4] == 'EPSG':
                proj.ImportFromEPSG(int(proj_string[5:]))
        gpkg.SetProjection(proj.ExportToWkt())
        gpkg.SetGeoTransform(geotransform)
        gpkg = None
    except Exception as e:
        sys.stderr.write(
            "ERROR: Cannot create SQLite GeoPackage '%s.gpkg'. "
            "Error message was: '%s'.\n" % (gpkg_name, e.message)
        )
        sys.exit(1)

def write_geotiff(fname, data, geo_transform, projection, data_type=gdal.GDT_Byte):
    """
    Create a GeoTIFF file with the given data.
    :param fname: Path to a directory with shapefiles
    :param data: Number of rows of the result
    :param geo_transform: Returned value of gdal.Dataset.GetGeoTransform (coefficients for
                          transforming between pixel/line (P,L) raster space, and projection
                          coordinates (Xp,Yp) space.
    :param projection: Projection definition string (Returned by gdal.Dataset.GetProjectionRef)
    """
    driver = gdal.GetDriverByName('GTiff')
    rows, cols = data.shape
    dataset = driver.Create(fname, cols, rows, 1, data_type)
    dataset.SetGeoTransform(geo_transform)
    dataset.SetProjection(projection)
    band = dataset.GetRasterBand(1)
    band.WriteArray(data)

    ct = gdal.ColorTable()
    for pixel_value in range(len(classes)+1):
        color_hex = COLORS[pixel_value]
        r = int(color_hex[1:3], 16)
        g = int(color_hex[3:5], 16)
        b = int(color_hex[5:7], 16)
        ct.SetColorEntry(pixel_value, (r, g, b, 255))
    band.SetColorTable(ct)

    metadata = {
        'TIFFTAG_COPYRIGHT': 'CC BY 4.0',
        'TIFFTAG_DOCUMENTNAME': 'classification',
        'TIFFTAG_IMAGEDESCRIPTION': 'Supervised classification.',
        'TIFFTAG_MAXSAMPLEVALUE': str(len(classes)),
        'TIFFTAG_MINSAMPLEVALUE': '0',
        'TIFFTAG_SOFTWARE': 'Python, GDAL, scikit-learn'
    }
    dataset.SetMetadata(metadata)

    dataset = None  # Close the file
    return

def array2raster(newRaster, RefRaster, array, noData, dataType):
    #data type conversion
    NP2GDAL_CONVERSION = { "uint8": 1, "int8": 1, "uint16": 2, "int16": 3, 
                          "uint32": 4, "int32": 5, "float32": 6, "float64": 7,
                          "complex64": 10, "complex128": 11,
                         }
    #get info from reference raster
    rfRaster = gdal.Open(RefRaster)
    geotransform = rfRaster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = array.shape[1]
    rows = array.shape[0]
    #create new raster
    outRaster = gdal.GetDriverByName('GTiff').Create(newRaster, cols, rows,1, NP2GDAL_CONVERSION[dataType])
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    #write array to band
    outband = outRaster.GetRasterBand(1)
    outband.SetNoDataValue(noData)
    outband.WriteArray(array)
    #define new raster projection
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(rfRaster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    #write raster
    outband.FlushCache()
    del rfRaster

def tif2png(inputFile,outputFile,verbose=True):
    '''Convert GeoTIFF to 8-bit RGB PNG'''
    src_ds = gdal.Open( inputFile ) 

    #Open output format driver, see gdal_translate --formats for list 
    format = "PNG"
    driver = gdal.GetDriverByName( format )

    ds = gdal.Open(inputFile)
    band1= ds.GetRasterBand(1)
    band2= ds.GetRasterBand(2)
    band3= ds.GetRasterBand(3)
    r = band1.ReadAsArray()
    g = band2.ReadAsArray()
    b = band3.ReadAsArray()
    #plt.figure(1)
    #plt.imshow(r)
    #plt.figure(2)
    #plt.imshow(g)
    #plt.figure(3)
    #plt.imshow(r)
    #plt.show()
    #Output to new format
    dst_ds = driver.CreateCopy( outputFile, src_ds, 1 )
    img = Image.open(outputFile)
    datas = img.getdata()
    newData = []
    y=0
    x=0
    for item in datas:
        rgb=(transform(r[x][y]),transform(g[x][y]),transform(b[x][y]))
        newData.append(rgb)
        if y<len(r[0])-1:
            y+=1
        else:
            if verbose==True:
                print str(x)+"\r",
            y=0
            x+=1
    img.putdata(newData)
    img.save(outputFile)

def create_raster(xsize, ysize, driver='MEM',tmpfile='', gt=None, srs_wkt=None, nodata=None,  init=None, datatype=gdal.GDT_Byte):
    # Create a memory raster to rasterize into.
    out_ds = gdal.GetDriverByName(driver).Create(tmpfile, xsize, ysize, 1 ,datatype)
    if init is not None:out_ds.GetRasterBand(1).Fill(init)
    if nodata is not None:out_ds.GetRasterBand(1).SetNoDataValue(nodata)
    if gt:out_ds.SetGeoTransform(gt)
    if srs_wkt:out_ds.SetProjection(srs_wkt)
    return out_ds

def get_ds(self):
        nl = self.ma_stack.shape[1]
        ns = self.ma_stack.shape[2]
        gdal_dtype = iolib.np2gdal_dtype(np.dtype(self.dtype))
        m_ds = gdal.GetDriverByName('MEM').Create('', ns, nl, 1, gdal_dtype)
        m_gt = [self.extent[0], self.res, 0, self.extent[3], 0, -self.res]
        m_ds.SetGeoTransform(m_gt)
        #this should already be WKT
        m_ds.SetProjection(self.proj)
        return m_ds

def write_datestack(self):
        #stat_list = ['dt_stack_ptp', 'dt_stack_mean', 'dt_stack_min', 'dt_stack_max', 'dt_stack_center']
        stat_list = ['dt_stack_ptp', 'dt_stack_min', 'dt_stack_max', 'dt_stack_center']
        if any([not hasattr(self, i) for i in stat_list]):
            #self.make_datestack()
            self.compute_dt_stats()
        print("Writing out datestack stats")
        #Create dummy ds - might want to use vrt here instead
        driver = gdal.GetDriverByName("MEM")
        ds = driver.Create('', self.dt_stack_ptp.shape[1], self.dt_stack_ptp.shape[0], 1, gdal.GDT_Float32)
        ds.SetGeoTransform(self.gt)
        ds.SetProjection(self.proj)
        #Write out with malib, should preserve ma type
        out_prefix = os.path.splitext(self.stack_fn)[0]
        iolib.writeGTiff(self.dt_stack_ptp, out_prefix+'_dt_ptp.tif', ds)
        self.dt_stack_ptp.set_fill_value(-9999)
        #iolib.writeGTiff(self.dt_stack_mean, out_prefix+'_dt_mean.tif', ds)
        #self.dt_stack_mean.set_fill_value(-9999)
        iolib.writeGTiff(self.dt_stack_min, out_prefix+'_dt_min.tif', ds)
        self.dt_stack_min.set_fill_value(-9999)
        iolib.writeGTiff(self.dt_stack_max, out_prefix+'_dt_max.tif', ds)
        self.dt_stack_max.set_fill_value(-9999)
        iolib.writeGTiff(self.dt_stack_center, out_prefix+'_dt_center.tif', ds)
        self.dt_stack_center.set_fill_value(-9999)

    #Note: want ot change the variable names min/max here
    #Might be better to save out as multiband GTiff here

def save(self, filename, driver=None, *args, **kwargs):

      if driver is None:
        ext = os.path.splitext(filename)[1][1:]
        if ext.lower() in ['tif', 'tiff']:
          driver = gdal.GetDriverByName('GTiff')
        else:
          raise Exception('Unkown extension. Can not determine driver')

      bands = self.array.shape[2] if len(self.array.shape)>2 else 1

      self.object = driver.Create(filename, self.array.shape[1], 
          self.array.shape[0], bands, GdalWriter.gdal_array_types[np.dtype(self.dtype)])

      if bands==1:
        self.object.GetRasterBand(1).WriteArray(self.array)
      else:
        for band in range(bands):
          self.object.GetRasterBand(band+1).WriteArray(self.array[:,:,band])

      #del self.object
      #Need to be deleted to actually save

      # dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )

      # srs = osr.SpatialReference()
      # srs.SetUTM( 11, 1 )
      # srs.SetWellKnownGeogCS( 'NAD27' )
      # dst_ds.SetProjection( srs.ExportToWkt() )



# @@@ Common feel functions @@@

def make_metatile(pile_directory, desired_value):
    metatile_filename = os.path.join(pile_directory,
                                     desired_value + '.tif')
    metatile_ds = gdal.GetDriverByName('GTiff').Create(
        metatile_filename, 4096, 4096, 1, gdal.GDT_Float32)

    # TODO: Try to add georeferencing...

    return metatile_filename, metatile_ds

def make_metatile(pile_directory):
    mosaic_filename = os.path.join(pile_directory,'mosaic.tif')
    mosaic_ds = gdal.GetDriverByName('GTiff').Create(
        mosaic_filename, 4096, 4096, 1, gdal.GDT_UInt16)

    # TODO: Try to add georeferencing...

    return mosaic_filename, mosaic_ds

def arr2Raster(data_name,Dsm_arr):
    [h,w] = Dsm_arr.shape
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(
            data_name, w, h, 1, gdal.GDT_Float32)
    dataset.GetRasterBand(1).WriteArray(Dsm_arr[:, :])
    dataset.FlushCache()

def open_vector(filename, driver=None):
    """Open vector file, return gdal.Dataset and OGR.Layer

        .. warning:: dataset and layer have to live in the same context,
            if dataset is deleted all layer references will get lost

        .. versionadded:: 0.12.0

    Parameters
    ----------
    filename : string
        vector file name
    driver : string
        gdal driver string

    Returns
    -------
    dataset : gdal.Dataset
        dataset
    layer : ogr.Layer
        layer
    """
    dataset = gdal.OpenEx(filename)

    if driver:
        gdal.GetDriverByName(driver)

    layer = dataset.GetLayer()

    return dataset, layer

def open_shape(filename, driver=None):
    """Open shapefile, return gdal.Dataset and OGR.Layer

        .. warning:: dataset and layer have to live in the same context,
            if dataset is deleted all layer references will get lost

        .. versionadded:: 0.6.0

    Parameters
    ----------
    filename : string
        shapefile name
    driver : string
        gdal driver string

    Returns
    -------
    dataset : gdal.Dataset
        dataset
    layer : ogr.Layer
        layer
    """

    if driver is None:
        driver = ogr.GetDriverByName('ESRI Shapefile')
    dataset = driver.Open(filename)
    if dataset is None:
        print('Could not open file')
        raise IOError
    layer = dataset.GetLayer()
    return dataset, layer

def open_raster(filename, driver=None):
    """Open raster file, return gdal.Dataset

        .. versionadded:: 0.6.0

    Parameters
    ----------
    filename : string
        raster file name
    driver : string
        gdal driver string

    Returns
    -------
    dataset : gdal.Dataset
        dataset
    """

    dataset = gdal.OpenEx(filename)

    if driver:
        gdal.GetDriverByName(driver)

    return dataset

def read_safnwc(filename):
    """Read MSG SAFNWC hdf5 file into a gdal georeferenced object

    Parameters
    ----------
    filename : string
        satellite file name

    Returns
    -------
    ds : gdal.DataSet
        with satellite data
    """

    root = gdal.Open(filename)
    ds1 = gdal.Open('HDF5:' + filename + '://CT')
    ds = gdal.GetDriverByName('MEM').CreateCopy('out', ds1, 0)

    # name = os.path.basename(filename)[7:11]
    try:
        proj = osr.SpatialReference()
        proj.ImportFromProj4(ds.GetMetadata()["PROJECTION"])
    except Exception:
        raise NameError("No metadata for satellite file %s" % filename)
    geotransform = root.GetMetadata()["GEOTRANSFORM_GDAL_TABLE"].split(",")
    geotransform[0] = root.GetMetadata()["XGEO_UP_LEFT"]
    geotransform[3] = root.GetMetadata()["YGEO_UP_LEFT"]
    ds.SetProjection(proj.ExportToWkt())
    ds.SetGeoTransform([float(x) for x in geotransform])
    return ds

def reprojectPanBand(src,match_ds,dst_filename):
    src_proj = src.GetProjection()
    src_geotrans = src.GetGeoTransform()

    match_proj = match_ds.GetProjection()
    match_geotrans = match_ds.GetGeoTransform()
    wide = match_ds.RasterXSize
    high = match_ds.RasterYSize

    dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Int16)
    dst.SetGeoTransform( match_geotrans )
    dst.SetProjection( match_proj)

    gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)
    return(gdal.Open(dst_filename,gdalconst.GA_ReadOnly))

def dump_raster(rast, rast_path, xoff=0, yoff=0, driver='GTiff', nodata=None):
    '''
    Creates a raster file from a given GDAL dataset (raster). Arguments:
        rast        A gdal.Dataset; does NOT accept NumPy array
        rast_path   The path of the output raster file
        xoff        Offset in the x-direction; should be provided when clipped
        yoff        Offset in the y-direction; should be provided when clipped
        driver      The name of the GDAL driver to use (determines file type)
        nodata      The NoData value; defaults to -9999.
    '''
    driver = gdal.GetDriverByName(driver)
    sink = driver.Create(rast_path, rast.RasterXSize, rast.RasterYSize,
        rast.RasterCount, rast.GetRasterBand(1).DataType)
    assert sink is not None, 'Cannot create dataset; there may be a problem with the output path you specified'
    sink.SetGeoTransform(rast.GetGeoTransform())
    sink.SetProjection(rast.GetProjection())

    for b in range(1, rast.RasterCount + 1):
        dat = rast.GetRasterBand(b).ReadAsArray()
        sink.GetRasterBand(b).WriteArray(dat)
        sink.GetRasterBand(b).SetStatistics(*map(np.float64,
            [dat.min(), dat.max(), dat.mean(), dat.std()]))

        if nodata is None:
            nodata = rast.GetRasterBand(b).GetNoDataValue()

            if nodata is None:
                nodata = -9999

        sink.GetRasterBand(b).SetNoDataValue(np.float64(nodata))

    sink.FlushCache()

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 reproject_tiff_custom(old_name, new_name, new_x_size, new_y_size,
                          new_geo_transform, new_projection, unit, mode):
    """
    Reprojects an tiff with custom tiff arguments. Can be used to warp tiff.
    If no projection is provided, fallback to old projection.
    Keyword arguments:
    old_name -- the name of the old tiff file
    new_name -- the name of the output tiff file
    new_x_size -- the number of new size in x dimension
    new_y_size -- the number of new size in y dimension
    new_geo_transform -- the new geo transform to apply
    new_projection -- the new projection to use
    unit -- the gdal unit in which the operation will be performed
    mode -- the gdal mode used for warping
    """
    mem_drv = gdal.GetDriverByName("MEM")

    old = gdal.Open(old_name)
    r = old.GetRasterBand(1)
    r.GetNoDataValue()

    # Adds 1 to keep the original zeros (reprojectImage maps NO_DATA to 0)
    old_array = old.ReadAsArray()
    mask = old_array == old.GetRasterBand(1).GetNoDataValue()
    old_array += 1
    old_array[mask] = 0
    temp = mem_drv.Create("temp", old.RasterXSize, old.RasterYSize, 1, unit)
    temp.SetGeoTransform(old.GetGeoTransform())
    temp.SetProjection(old.GetProjection())
    temp.GetRasterBand(1).WriteArray(old_array)

    new = mem_drv.Create(new_name, new_x_size, new_y_size, 1, unit)

    new.SetGeoTransform(new_geo_transform)
    if new_projection is None:
        new.SetProjection(old.GetProjection())
    else:
        new.SetProjection(new_projection)

    res = gdal.ReprojectImage(temp, new, old.GetProjection(), new_projection,
                              mode)

    assert res == 0, 'Error in ReprojectImage'
    arr = new.ReadAsArray()
    mask = arr != 0
    arr -= 1
    arr[~mask] = 0
    new = None
    temp = None

    return arr, mask

def makeDensityRaster(speciesOcc, inVector, pixelSize, outRas, noData):
    srcVector = ogr.Open(inVector)
    srcLayer = srcVector.GetLayer()
    srs = srcLayer.GetSpatialRef()
    # if the layer is not wgs84
    if srs.GetAttrValue("AUTHORITY", 1) != '4326':
        print('Layer projection should be WGS84!')
        return

    xMin, xMax, yMin, yMax = srcLayer.GetExtent()

    # Create the destination data source
    xRes = int((xMax - xMin) / pixelSize)
    yRes = int((yMax - yMin) / pixelSize)
    targetRas = gdal.GetDriverByName('GTiff').Create(outRas, xRes, yRes, 1, 6) # 6 == float
    targetRas.SetGeoTransform((xMin, pixelSize, 0, yMax, 0, -pixelSize))
    band = targetRas.GetRasterBand(1)
    band.SetNoDataValue(noData)

    # Rasterize clips the raster band
    gdal.RasterizeLayer(targetRas, [1], srcLayer, None, None, [0], ['ALL_TOUCHED=TRUE'])
    #os.remove(outRas)
    g = band.ReadAsArray()

    for point in speciesOcc.iterrows():
        xi = int((point[1]['x'] - xMin) / pixelSize)
        yi = int((point[1]['y'] - yMax) / -pixelSize)

        try:
            if g[yi,xi] != noData:
                g[yi,xi] += 1
            else:
                print('point ({}, {}) out of bounds'.format(point[1]['x'], point[1]['y']))
        except:
            print('point ({}, {}) out of bounds'.format(point[1]['x'], point[1]['y']))
            pass


    band.WriteArray(g)
    targetRasSRS = osr.SpatialReference()
    targetRasSRS.ImportFromWkt(srs.ExportToWkt())
    targetRas.SetProjection(targetRasSRS.ExportToWkt())
    band.FlushCache()

def createRaster(outRas, extCells, pixelSize, cellValues = 'lat', dataType = "float32", noData = -9999, inVector = None):
    NP2GDAL_CONVERSION = { "uint8": 1, "uint16": 2, "int16": 3, 
                          "uint32": 4, "int32": 5, "float32": 6, "float64": 7,
                          "complex64": 10, "complex128": 11,
                         }
    if os.path.exists(outRas):
        print('Raster file already excists!')
        return
    # Create the destination data source        
    xRes = int(numpy.ceil((extCells[2] - extCells[0]) / pixelSize))
    yRes = int(numpy.ceil((extCells[3] - extCells[1]) / pixelSize))

    targetRas = gdal.GetDriverByName('GTiff').Create(outRas, xRes, yRes, 1, NP2GDAL_CONVERSION[dataType])
    targetRas.SetGeoTransform((extCells[0], pixelSize, 0, extCells[3], 0, -pixelSize))
    band = targetRas.GetRasterBand(1)
    band.SetNoDataValue(noData)

    if inVector != None:
        srcVector = ogr.Open(inVector)
        srcLayer = srcVector.GetLayer()
        srs = srcLayer.GetSpatialRef()
        # if the layer is not wgs84
        if srs.GetAttrValue("AUTHORITY", 1) != '4326':
            print('Layer projection should be WGS84!')
            os.remove(outRas)
            return

        # Rasterize clips the raster band
        gdal.RasterizeLayer(targetRas, [1], srcLayer, None, None, [0], ['ALL_TOUCHED=TRUE'])

        g = band.ReadAsArray()

    else:
        g = numpy.zeros((yRes,xRes), eval('numpy.{}'.format(dataType)))

    #populate matrix with random numbers
    for i in range(yRes):
        for j in range(xRes):
            if g[i,j] != noData:
                if cellValues == 'lat':
                    g[i,j] = i
                elif cellValues == 'lon':
                    g[i,j] = j
                elif cellValues == 'random':
                    g[i,j] = numpy.random.randint(50)

    band.WriteArray(g)
    targetRasSRS = osr.SpatialReference()
    targetRasSRS.ImportFromEPSG(4326)
    targetRas.SetProjection(targetRasSRS.ExportToWkt())
    band.FlushCache()

#function to filter raster cells based on the coverage by some vector features

def geom2shp(geom, out_fn, fields=False):
    """Write out a new shapefile for input geometry
    """
    from pygeotools.lib import timelib
    driverName = "ESRI Shapefile"
    drv = ogr.GetDriverByName(driverName)
    if os.path.exists(out_fn):
        drv.DeleteDataSource(out_fn)
    out_ds = drv.CreateDataSource(out_fn)
    out_lyrname = os.path.splitext(os.path.split(out_fn)[1])[0]
    geom_srs = geom.GetSpatialReference()
    geom_type = geom.GetGeometryType()
    out_lyr = out_ds.CreateLayer(out_lyrname, geom_srs, geom_type)
    if fields:
        field_defn = ogr.FieldDefn("name", ogr.OFTString)
        field_defn.SetWidth(128)
        out_lyr.CreateField(field_defn)
        field_defn = ogr.FieldDefn("path", ogr.OFTString)
        field_defn.SetWidth(254)
        out_lyr.CreateField(field_defn)
        #field_defn = ogr.FieldDefn("date", ogr.OFTString)
        #This allows sorting by date
        field_defn = ogr.FieldDefn("date", ogr.OFTInteger)
        field_defn.SetWidth(32)
        out_lyr.CreateField(field_defn)
        field_defn = ogr.FieldDefn("decyear", ogr.OFTReal)
        field_defn.SetPrecision(8)
        field_defn.SetWidth(64)
        out_lyr.CreateField(field_defn)
    out_feat = ogr.Feature(out_lyr.GetLayerDefn())
    out_feat.SetGeometry(geom)
    if fields:
        #Hack to force output extesion to tif, since out_fn is shp
        out_path = os.path.splitext(out_fn)[0] + '.tif'
        out_feat.SetField("name", os.path.split(out_path)[-1])
        out_feat.SetField("path", out_path)
        #Try to extract a date from input raster fn
        out_feat_date = timelib.fn_getdatetime(out_fn)
        if out_feat_date is not None:
            datestamp = int(out_feat_date.strftime('%Y%m%d'))
            #out_feat_date = int(out_feat_date.strftime('%Y%m%d%H%M'))
            out_feat.SetField("date", datestamp)
            decyear = timelib.dt2decyear(out_feat_date)
            out_feat.SetField("decyear", decyear)
    out_lyr.CreateFeature(out_feat)
    out_ds = None
    #return status?

def from_array(cls, raster, geo_transform, proj4, 
        gdal_dtype=gdal.GDT_Float32, nodata=None):
        """ Create a georaster object from numpy array and georeferencing information.

        :param raster: 2-D NumPy array of raster to load
        :type raster: np.array
        :param geo_transform: a Geographic Transform tuple of the form \
        (top left x, w-e cell size, 0, top left y, 0, n-s cell size (-ve))
        :type geo_transform: tuple
        :param proj4: a proj4 string representing the raster projection
        :type proj4: str
        :param gdal_dtype: a GDAL data type (default gdal.GDT_Float32)
        :type gdal_dtype: int
        :param nodata: None or the nodata value for this array
        :type nodata: None, int, float, str

        :returns: GeoRaster object
        :rtype: GeoRaster

         """

        if len(raster.shape) > 2:
            nbands = raster.shape[2]
        else:
            nbands = 1

        # Create a GDAL memory raster to hold the input array
        mem_drv = gdal.GetDriverByName('MEM')
        source_ds = mem_drv.Create('', raster.shape[1], raster.shape[0],
                         nbands, gdal_dtype)

        # Set geo-referencing
        source_ds.SetGeoTransform(geo_transform)
        srs = osr.SpatialReference()
        srs.ImportFromProj4(proj4)
        source_ds.SetProjection(srs.ExportToWkt())

        # Write input array to the GDAL memory raster
        for b in range(0,nbands):
            if nbands > 1:
                r = raster[:,:,b]
            else:
                r = raster
            source_ds.GetRasterBand(b+1).WriteArray(r)
            if nodata != None:
                source_ds.GetRasterBand(b+1).SetNoDataValue(nodata)

        # Return a georaster instance instantiated by the GDAL raster
        return cls(source_ds)

def creaLayer(self):
        layer = self.iface.activeLayer()

        provider = layer.dataProvider()

        path = provider.dataSourceUri()

        (raiz, filename) = os.path.split(path)

        dataset = gdal.Open(path)

        # Get projection
        prj = dataset.GetProjection()

        # setting band
        number_band = 1

        band = dataset.GetRasterBand(number_band)

        # Get raster metadata
        geotransform = dataset.GetGeoTransform()

        # Set name of output raster
        output_file = "C:\\Users\\caligola\\Desktop\\s\\raster_output.tif"

        # Create gtif file with rows and columns from parent raster
        driver = gdal.GetDriverByName("GTiff")

        raster = self.changeRasterValues(band)

        dst_ds = driver.Create(output_file,
                               band.XSize,
                               band.YSize,
                               number_band,
                               band.DataType)

        # writting output raster
        dst_ds.GetRasterBand(number_band).WriteArray(raster)

        # setting extension of output raster
        # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
        dst_ds.SetGeoTransform(geotransform)

        # setting spatial reference of output raster
        srs = osr.SpatialReference(wkt=prj)
        dst_ds.SetProjection(srs.ExportToWkt())

def reproject_dataset_example(dataset, dataset_example, method=1):

    # open dataset that must be transformed 
    try:
        if dataset.split('.')[-1] == 'tif':
            g = gdal.Open(dataset)                                  
        else:
            g = dataset    
    except:
            g = dataset            
    epsg_from = Get_epsg(g)    

    # open dataset that is used for transforming the dataset
    try:
        if dataset_example.split('.')[-1] == 'tif':
            gland = gdal.Open(dataset_example)                                  
        else:
            gland = dataset_example      
    except:
            gland = dataset_example              
    epsg_to = Get_epsg(gland)   

    # Set the EPSG codes
    osng = osr.SpatialReference()
    osng.ImportFromEPSG(epsg_to)
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(epsg_from)

    # Get shape and geo transform from example              
    geo_land = gland.GetGeoTransform()          
    col=gland.RasterXSize
    rows=gland.RasterYSize

    # Create new raster         
    mem_drv = gdal.GetDriverByName('MEM')
    dest1 = mem_drv.Create('', col, rows, 1, gdal.GDT_Float32)
    dest1.SetGeoTransform(geo_land)
    dest1.SetProjection(osng.ExportToWkt())

    # Perform the projection/resampling
    if method is 1:             
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_NearestNeighbour)
    if method is 2:             
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Bilinear)
    if method is 3:             
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Lanczos)
    if method is 4:             
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Average)
    return(dest1)