Python osgeo.gdal.GetDriverByName() Examples

def make_raster(in_ds, fn, data, data_type, nodata=None):
    """Create a one-band GeoTIFF.

    in_ds     - datasource to copy projection and geotransform from
    fn        - path to the file to create
    data      - NumPy array containing data to write
    data_type - output data type
    nodata    - optional NoData value
    """
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(
        fn, in_ds.RasterXSize, in_ds.RasterYSize, 1, data_type)
    out_ds.SetProjection(in_ds.GetProjection())
    out_ds.SetGeoTransform(in_ds.GetGeoTransform())
    out_band = out_ds.GetRasterBand(1)
    if nodata is not None:
        out_band.SetNoDataValue(nodata)
    out_band.WriteArray(data)
    out_band.FlushCache()
    out_band.ComputeStatistics(False)
    return out_ds 

def createRaster(self, srcraster, destfile, pixelType):
        """
        Creates a raster file
        srcraster: source raster
        destfile: destination file
        pixelType: raster pixel type (e.g byte)
        """
        cols = srcraster.RasterXSize
        rows = srcraster.RasterYSize

        (xOrigin, yOrigin, pixelWidth, pixelHeight) = self.getGeoreferenceInfo(srcraster)

        targetSR = self.getCRS(srcraster)

        driver = gdal.GetDriverByName('GTiff')

        outRaster = driver.Create(destfile, cols, rows, 3, pixelType)
        outRaster.SetGeoTransform((xOrigin, pixelWidth, 0, yOrigin, 0, pixelHeight))
        outRaster.SetProjection(targetSR.ExportToWkt())
        
        return outRaster 

def createRaster(self, srcraster, destfile, pixelType):
        '''
        Creates a raster file
        srcraster: source raster
        destfile: destination file
        pixelType: raster pixel type (e.g byte)
        '''
        cols = srcraster.RasterXSize
        rows = srcraster.RasterYSize

        (xOrigin, yOrigin, pixelWidth, pixelHeight) = self.getGeoreferenceInfo(srcraster)

        targetSR = self.getCRS(srcraster)

        driver = gdal.GetDriverByName('GTiff')

        outRaster = driver.Create(destfile, cols, rows, 3, pixelType)
        outRaster.SetGeoTransform((xOrigin, pixelWidth, 0, yOrigin, 0, pixelHeight))
        outRaster.SetProjection(targetSR.ExportToWkt())
        
        return outRaster 

def CreateMultiBandGeoTiff(OutPath, Array):
    '''
    Author: Jake Shermeyer
    Array has shape:
        Channels, Y, X?
    '''
    driver = gdal.GetDriverByName('GTiff')
    DataSet = driver.Create(OutPath, Array.shape[2], Array.shape[1],
                            Array.shape[0], gdal.GDT_Byte,
                            ['COMPRESS=LZW'])
    for i, image in enumerate(Array, 1):
        DataSet.GetRasterBand(i).WriteArray(image)
    del DataSet

    return OutPath


############################################################################### 

def write_data(self, data, fileName, colorTable):
		fileName 	= os.path.join(self.outpath, fileName)
		driver 		= gdal.GetDriverByName( "GTiff" )
		dst_ds 		= driver.Create( fileName, self.RasterXSize, self.RasterYSize, 1, gdal.GDT_Byte, [ 'COMPRESS=DEFLATE' ] )
		band 		= dst_ds.GetRasterBand(1)

		if self.geotransform:
			dst_ds.SetGeoTransform( self.geotransform )
			
		if self.projection:
			dst_ds.SetProjection( self.projection )

		if colorTable:
			print "Add colortable"
			band.SetRasterColorTable(self.ct)

		band.WriteArray(data, 0, 0)
			
		if verbose:
			print "Written", fileName

		ds 		= None 

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 write_geotiff(tiff_path, data, upper_left_lon, upper_left_lat, step, AREA_OR_POINT='Point'):
    """
    写入geotiff。默认pixel is point，默认WGS84
    """
    rows, columns = data.shape
    bands = 1
    pixel_width = step
    pixel_height = -step
    half_step = step / 2
    upper_left_x = upper_left_lon - half_step
    upper_left_y = upper_left_lat + half_step
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(tiff_path, columns, rows, bands, gdal.GDT_Float32)
    dataset.SetMetadataItem('AREA_OR_POINT', AREA_OR_POINT)
    dataset.SetGeoTransform([upper_left_x, pixel_width, 0,
                             upper_left_y, 0, pixel_height])
    # 获取地理坐标系统信息，用于选取需要的地理坐标系统
    sr = osr.SpatialReference()
    sr.SetWellKnownGeogCS('WGS84')
    dataset.SetProjection(sr.ExportToWkt())  # 给新建图层赋予投影信息
    dataset.GetRasterBand(1).WriteArray(data)
    dataset.FlushCache()  # 将数据写入硬盘
    dataset = None 

def to_tiff(self, filename):
        '''
        geo.to_tiff(filename)

        Saves GeoRaster as geotiff filename.tif with type datatype

        If GeoRaster does not have datatype, then it tries to assign a type.
        You can assign the type yourself by setting
         geo.datatype = 'gdal.GDT_'+type
        '''
        if self.datatype is None:
            self.datatype = gdal_array.NumericTypeCodeToGDALTypeCode(self.raster.data.dtype)
            if self.datatype is None:
                if self.raster.data.dtype.name.find('int') !=- 1:
                    self.raster = self.raster.astype(np.int32)
                    self.datatype = gdal_array.NumericTypeCodeToGDALTypeCode(self.raster.data.dtype)
                else:
                    self.raster = self.raster.astype(np.float64)
                    self.datatype = gdal_array.NumericTypeCodeToGDALTypeCode(self.raster.data.dtype)
        self.raster.data[self.raster.mask] = self.nodata_value
        create_geotiff(filename, self.raster, gdal.GetDriverByName('GTiff'), self.nodata_value,
                       self.shape[1], self.shape[0], self.geot, self.projection, self.datatype) 

def make_tif(path, tloffset=(0, 0), reso=(0.25, -0.25), rsize=(20, 10),
             proj=SRS[0]['wkt'], channel_count=1, dtype=gdal.GDT_Float32):
    """Create a tiff files and return info about it"""
    tl = ROOT_TL + tloffset
    reso = np.asarray(reso)
    fp = buzz.Footprint(tl=tl, rsize=rsize, size=np.abs(reso * rsize))
    x, y = fp.meshgrid_spatial
    x = np.abs(x) - abs(ROOT_TL[0])
    y = abs(ROOT_TL[1]) - np.abs(y)
    x *= 15
    y *= 15
    a = x / 2 + y / 2
    a = np.around(a).astype('float32')
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(path, rsize[0], rsize[1], channel_count, dtype)
    dataset.SetGeoTransform(fp.gt)
    dataset.SetProjection(proj)
    for i in range(channel_count):
        dataset.GetRasterBand(i + 1).WriteArray(a)
        dataset.GetRasterBand(i + 1).SetNoDataValue(-32000.)
    dataset.FlushCache()
    return path, fp, a 

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

    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 main():
    args = parse_arguments()
    src_ds = gdal.Open(args.src_file[0])

    driver = gdal.GetDriverByName("GTiff")

    dst_ds = driver.CreateCopy(args.out_file[0], src_ds, 0)

    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(float(args.position[0]), float(args.position[1]))
    point.Transform(transform)

    gt = [point.GetX(), float(args.pixelsize[0]), 0,
          point.GetY(), 0, -float(args.pixelsize[1])]

    dst_ds.SetGeoTransform(gt)

    to_sys = osr.SpatialReference()
    to_sys.ImportFromEPSG(to_srs)

    dest_wkt = to_sys.ExportToWkt()

    dst_ds.SetProjection(dest_wkt) 

def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)
    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# convert images in a selected folder to shapefiles 

def toGDAL(self):
		'''Get GDAL memory driver dataset'''
		w, h = self.size
		n = self.nbBands
		dtype = str(self.dtype)
		if dtype == 'uint8': dtype = 'byte'
		dtype = gdal.GetDataTypeByName(dtype)
		mem = gdal.GetDriverByName('MEM').Create('', w, h, n, dtype)
		#writearray is available only at band level
		if self.isOneBand:
			mem.GetRasterBand(1).WriteArray(self.data)
		else:
			for bandIdx in range(n):
				bandArray = self.data[:,:,bandIdx]
				mem.GetRasterBand(bandIdx+1).WriteArray(bandArray)
		#write georef
		if self.isGeoref:
			mem.SetGeoTransform(self.georef.toGDAL())
			if self.georef.crs is not None:
				mem.SetProjection(self.georef.crs.getOgrSpatialRef().ExportToWkt())
		return mem 

def writeRaster(arr, out_path, template):
    no_data = 0
    # First of all, gather some information from the template file
    data = gdal.Open(template)
    [cols, rows] = arr.shape
    trans = data.GetGeoTransform()
    proj = data.GetProjection()
    # nodatav = 0 #data.GetNoDataValue()
    # Create the file, using the information from the template file
    outdriver = gdal.GetDriverByName("GTiff")
    # http://www.gdal.org/gdal_8h.html
    # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6,
    outdata   = outdriver.Create(str(out_path), rows, cols, 1, gdal.GDT_UInt32)
    # Write the array to the file, which is the original array in this example
    outdata.GetRasterBand(1).WriteArray(arr)
    # Set a no data value if required
    outdata.GetRasterBand(1).SetNoDataValue(no_data)
    # Georeference the image
    outdata.SetGeoTransform(trans)
    # Write projection information
    outdata.SetProjection(proj)
    return arr


# raster to vector 

def fillNodata(self):
		#if not self.noData in self.data:
		if not np.ma.is_masked(self.data):
			#do not process it if its not necessary
			return
		if self.IFACE == 'GDAL':
			# gdal.FillNodata need a band object to apply on
			# so we create a memory datasource (1 band, float)
			height, width = self.data.shape
			ds = gdal.GetDriverByName('MEM').Create('', width, height, 1, gdal.GetDataTypeByName('float32'))
			b = ds.GetRasterBand(1)
			b.SetNoDataValue(self.noData)
			self.data =  np.ma.filled(self.data, self.noData)# Fill mask with nodata value
			b.WriteArray(self.data)
			gdal.FillNodata(targetBand=b, maskBand=None, maxSearchDist=max(self.size.xy), smoothingIterations=0)
			self.data = b.ReadAsArray()
			ds, b = None, None
		else: #Call the inpainting function
			# Cast to float
			self.cast2float()
			# Fill mask with NaN (warning NaN is a special value for float arrays only)
			self.data =  np.ma.filled(self.data, np.NaN)
			# Inpainting
			self.data = replace_nans(self.data, max_iter=5, tolerance=0.5, kernel_size=2, method='localmean') 

def save(self, path):
		'''
		save the numpy array to a new image file
		output format is defined by path extension
		'''

		imgFormat = path[-3:]

		if self.IFACE == 'PIL':
			self.toPIL().save(path)
		elif self.IFACE == 'IMGIO':
			if imgFormat == 'jpg' and self.hasAlpha:
				self.removeAlpha()
			imageio.imwrite(path, self.data)#float32 support ok
		elif self.IFACE == 'GDAL':
			if imgFormat == 'png':
				driver = 'PNG'
			elif imgFormat == 'jpg':
				driver = 'JPEG'
			elif imgFormat == 'tif':
				driver = 'Gtiff'
			else:
				raise ValueError('Cannot write to '+ imgFormat + ' image format')
			#Some format like jpg or png has no create method implemented
			#because we can't write data at random with these formats
			#so we must use an intermediate memory driver, write data to it
			#and then write the output file with the createcopy method
			mem = self.toGDAL()
			out = gdal.GetDriverByName(driver).CreateCopy(path, mem)
			mem = out = None

		if self.isGeoref:
			self.georef.toWorldFile(os.path.splitext(path)[0] + '.wld') 

def writeGTiff(dstImage, outputPath, projection, coordinates, dtype = None):
    """
    @brief Writes a GeoTiff file created from an existing coordinate-system
    @Note Coordinates: [Top-Left X, W-E Resolution, 0, Top Left Y, 0, N-S Resolution]
    """
    driver = gdal.GetDriverByName('GTiff')
    #Add dimension for a single band-image
    if(dstImage.ndim is 2):
        dstImage = dstImage[..., np.newaxis]  
    #Set output dtype if not specified
    if(dtype == None):
        dtype = gdal_array.NumericTypeCodeToGDALTypeCode(dstImage.dtype)
    dataset = driver.Create(outputPath, dstImage.shape[1], dstImage.shape[0] , dstImage.shape[2], dtype)
    if(dataset == None):
        raise OSError ("GDAL Could not create {0}".format(outputPath))
    if(not coordinates or len(coordinates) != 6):
        raise ValueError("Geotransform empty or unknown: {0}".format(coordinates))
    dataset.SetGeoTransform(coordinates)  
    if(not projection):
        raise ValueError("Projection empty or unknown: {0}".format(projection))
    dataset.SetProjection(projection)
    for index in range(dstImage.shape[2]):
        dataset.GetRasterBand(index+1).WriteArray(dstImage[:,:,index])
    dataset.FlushCache()  # Write to disk.
    dataset = None
    driver = None
    return 0 

def writer(self, array, *args, **kwargs):
        if hasattr(pyrat, "app"):  # select correct GDAL driver
            drv_idx = self.dnames.index(self.file[1])  # GUI mode: long name
            driver = gdal.GetDriverByName(self.driver[drv_idx])
        else:
            driver = gdal.GetDriverByName(self.file[1])  # CLI mode: short name

        if driver is None:
            logging.error("Unknown GDAL driver: " + self.file[1])
            return

        meta = kwargs["meta"]

        shp = array.shape  # reshape to (nchannel, ny, nx)
        array = array.reshape((int(np.prod(shp[0:-2])),) + shp[-2:])
        shp = array.shape
        gdaltype = NP2GDAL_CONVERSION[array.dtype.name]

        ds = driver.Create(self.file[0], shp[2], shp[1], shp[0], gdaltype)
        for k, arr in enumerate(array):
            ds.GetRasterBand(k + 1).WriteArray(arr)

        gmeta, cmeta = dict2gdalmeta(meta, shp[0])
        ds.SetMetadata(gmeta, "")
        for k in range(ds.RasterCount):
            ds.GetRasterBand(k + 1).SetMetadata(cmeta[k], "")

        # ------------------- Comment: No handling of projections yet
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)
        ds.SetProjection(srs.ExportToWkt())
        # -------------------

        ds = None  # correct according to GDAL manual!!?? 

def TestLand(lon, lat):
    latlon = osr.SpatialReference()
    latlon.ImportFromEPSG(4326)

    # create a point

    pt = ogr.Geometry(ogr.wkbPoint)
    pt.SetPoint_2D(0, lon, lat)

    # read shapefile
    shapefile = "land_polygons_osm/simplified_land_polygons.shp"
    driver = ogr.GetDriverByName("ESRI Shapefile")
    dataSource = driver.Open(shapefile, 0)
    layer = dataSource.GetLayer()
    targetProj = layer.GetSpatialRef()
    land = False

    # conversion to shapefile projection
    transform = osr.CoordinateTransformation(latlon, targetProj)
    pt.Transform(transform)

    # search point in layers
    for feature in layer:
        geom = feature.GetGeometryRef()
        if geom.Contains(pt):
            land = True
            break

    return land


##################################### Lecture de fichier de parametres "Mot_clé=Valeur" 

def toBLOB(self, ext='PNG'):
		'''Get bytes raw data'''
		if ext == 'JPG':
			ext = 'JPEG'

		if self.IFACE == 'PIL':
			b = io.BytesIO()
			img = Image.fromarray(self.data)
			img.save(b, format=ext)
			data = b.getvalue() #convert bytesio to bytes

		elif self.IFACE == 'IMGIO':
			if ext == 'JPEG' and self.hasAlpha:
				self.removeAlpha()
			data = imageio.imwrite(imageio.RETURN_BYTES, self.data, format=ext)

		elif self.IFACE == 'GDAL':
			mem = self.toGDAL()
			#build a random name to make the function thread safe
			name = ''.join(random.choice('abcdefghijklmnopqrstuvwxyz') for i in range(5))
			vsiname = '/vsimem/' + name + '.png'
			out = gdal.GetDriverByName(ext).CreateCopy(vsiname, mem)
			# Read /vsimem/output.png
			f = gdal.VSIFOpenL(vsiname, 'rb')
			gdal.VSIFSeekL(f, 0, 2) # seek to end
			size = gdal.VSIFTellL(f)
			gdal.VSIFSeekL(f, 0, 0) # seek to beginning
			data = gdal.VSIFReadL(1, size, f)
			gdal.VSIFCloseL(f)
			# Cleanup
			gdal.Unlink(vsiname)
			mem = None

		return data 

def _crop_resample_setup(extents, input_tif, new_res, output_file,
                         dst_driver_type='GTiff', out_bands=2):
    """
    Convenience function for crop/resample setup
    """
    # Source
    src_ds = gdal.Open(input_tif, gdalconst.GA_ReadOnly)
    src_proj = src_ds.GetProjection()

    # source metadata to be copied into the output
    meta_data = src_ds.GetMetadata()

    # get the image extents
    min_x, min_y, max_x, max_y = extents
    geo_transform = src_ds.GetGeoTransform()  # tuple of 6 numbers

    # Create a new geotransform for the image
    gt2 = list(geo_transform)
    gt2[0] = min_x
    gt2[3] = max_y
    # We want a section of source that matches this:
    resampled_proj = src_proj
    if new_res[0]:  # if new_res is not None, it can't be zero either
        resampled_geotrans = gt2[:1] + [new_res[0]] + gt2[2:-1] + [new_res[1]]
    else:
        resampled_geotrans = gt2

    px_height, px_width = _gdalwarp_width_and_height(max_x, max_y, min_x,
                                                     min_y, resampled_geotrans)

    # Output / destination
    dst = gdal.GetDriverByName(dst_driver_type).Create(
        output_file, px_width, px_height, out_bands, gdalconst.GDT_Float32)
    dst.SetGeoTransform(resampled_geotrans)
    dst.SetProjection(resampled_proj)

    for k, v in meta_data.items():
        dst.SetMetadataItem(k, v)

    return dst, resampled_proj, src_ds, src_proj 

def _writeRaster(self, data, filename):
        """Writes GeoTIFF raster temporarily so that it can be imported into the database."""
        nrows, ncols = data.shape
        driver = gdal.GetDriverByName("GTiff")
        ods = driver.Create(filename, ncols, nrows, 1, gdal.GDT_Float32)
        ods.SetGeoTransform([min(self.lon) - self.res / 2.0, self.res,
                             0, max(self.lat) + self.res / 2.0, 0, -self.res])
        srs = osr.SpatialReference()
        srs.SetWellKnownGeogCS("WGS84")
        ods.SetProjection(srs.ExportToWkt())
        ods.GetRasterBand(1).WriteArray(data)
        ods.GetRasterBand(1).SetNoDataValue(self.nodata)
        ods = None 

def rgb2gray( self, red, green, blue, RasterXSize, RasterYSize ):		
		lum			= 0.2126 * red + 0.7152 * green + 0.0722 * blue
		
		#driver 		= gdal.GetDriverByName( "GTiff" )
		#dst_ds 		= driver.Create( self.output_file2, RasterXSize, RasterYSize, 1, gdal.GDT_Byte, [ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )
		#band 		= dst_ds.GetRasterBand(1)
		#band.WriteArray(lum, 0, 0)
		#band.SetNoDataValue(0)
		#dst_ds 		= None
		#print "Written", self.output_file2
		return lum 

def pts2raster(shapefile,RASTER_PATH,cellSize,field_name=False):
    from osgeo import gdal, ogr

    # Define pixel_size and NoData value of new raster
    pixel_size = cellSize
    NoData_value = -9999
    
    # Filename of input OGR file
    vector_ptsShp_fn = shapefile
    
    # Filename of the raster Tiff that will be created
    raster_ptsShp_fn = RASTER_PATH
    
    # Open the data source and read in the extent
    source_ds = ogr.Open(vector_ptsShp_fn)
    source_layer = source_ds.GetLayer()
    x_min, x_max, y_min, y_max = source_layer.GetExtent()
    
    # Create the destination data source
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_ptsShp_fn, x_res, y_res, 1, gdal.GDT_Int32 )
    target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(NoData_value)
    
    # Rasterize
    # gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])
    # Rasterize
    if field_name:
        gdal.RasterizeLayer(target_ds,[1], source_layer,options=["ATTRIBUTE={0}".format(field_name)])
        # print("write:",field_name)
    else:
        gdal.RasterizeLayer(target_ds,[1], source_layer,burn_values=[-1])   
    return gdal.Open(RASTER_PATH).ReadAsArray() 

#批量计算 

def write_data(self, data, fileName, type=gdal.GDT_Float32, nbands=1, ct=0):
		fileName 	= os.path.join(self.outpath, fileName)
		
		if self.verbose:
			print "write_data", fileName
		
		driver 		= gdal.GetDriverByName( "GTiff" )
		dst_ds 		= driver.Create( fileName, self.RasterXSize, self.RasterYSize, nbands, type, [ 'INTERLEAVE=PIXEL', 'COMPRESS=DEFLATE' ] )
		band 		= dst_ds.GetRasterBand(1)

		band.WriteArray(data, 0, 0)
		band.SetNoDataValue(0)

		if ct :
			if self.verbose:
				print "write ct"
			
			ct = self.generate_color_table()
			band.SetRasterColorTable(ct)
			
		if self.verbose:
			print "write geotransform and projection"
		dst_ds.SetGeoTransform( self.geotransform )
		dst_ds.SetProjection( self.projection )
			
		if self.verbose:
			print "Written", fileName

		dst_ds 		= None 

def CreateTopojsonFile(self, fileName, src_ds, projection, geotransform, ct, data, pres, xorg, ymax, water ):
	
		driver 				= gdal.GetDriverByName( "GTiff" )
		dst_ds_dataset		= driver.Create( fileName, src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Byte, [ 'COMPRESS=DEFLATE' ] )
		
		dst_ds_dataset.SetGeoTransform( geotransform )
		dst_ds_dataset.SetProjection( projection )

		o_band				= dst_ds_dataset.GetRasterBand(1)
	
		o_band.SetRasterColorTable(ct)
		o_band.WriteArray(data, 0, 0)

		dst_ds_dataset = None
		print "Created", fileName

		cmd = "gdal_translate -q -of PNM -expand gray " + fileName + " "+fileName+".bmp"
		self.execute(cmd)

		# -i  		invert before processing
		# -t 2  	suppress speckles of up to this many pixels. 
		# -a 1.5  	set the corner threshold parameter
		# -z black  specify how to resolve ambiguities in path decomposition. Must be one of black, white, right, left, minority, majority, or random. Default is minority
		# -x 		scaling factor
		# -L		left margin
		# -B		bottom margin

		cmd = str.format("potrace -i -z black -a 1.5 -t 3 -b geojson -o {0} {1} -x {2} -L {3} -B {4} ", fileName+".geojson", fileName+".bmp", pres, xorg, ymax ); 
		self.execute(cmd)
	
		cmd = str.format("topojson -o {0} --simplify-proportion 0.75 -p water={1} -- water={2}", fileName+".topojson", water, fileName+".geojson" ); 
		self.execute(cmd)
	
		# convert it back to json
		cmd = "topojson-geojson --precision 4 -o %s %s" % ( self.geojsonDir, fileName+".topojson" )
		self.execute(cmd)
	
		# rename file
		output_file = "water_level_%d.geojson" % water
		cmd = "mv %s %s" % (os.path.join(self.geojsonDir,"water.json"), os.path.join(self.geojsonDir, output_file))
		self.execute(cmd) 

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)

    try:
        proj = osr.SpatialReference()
        proj.ImportFromProj4(ds.GetMetadata()["PROJECTION"])
    except KeyError:
        raise KeyError(
            "WRADLIB: Projection is missing for satellite " "file {}".format(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 open_vector(filename, driver=None, layer=0):
    """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

    Parameters
    ----------
    filename : string
        vector file name
    driver : string
        gdal driver string
    layer : int or string

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

    if driver:
        gdal.GetDriverByName(driver)

    layer = dataset.GetLayer(layer)

    return dataset, layer 

def update_msi(input_file_name, output_file_name):
    img = tifffile.imread(input_file_name)
    rows, cols, bands = img.shape
    driver = gdal.GetDriverByName("GTiff")
    output_data = driver.Create(output_file_name, rows, cols, bands, gdal.GDT_UInt16)
    for band in range(0, bands):
        output_band = output_data.GetRasterBand(band + 1)
        output_band.WriteArray(img[:, :, band])
    output_data.FlushCache()
    output_data = None


# main 

def write(arr, filename):
        """
        Write image array to a file with the given filename.

        Args:
            arr (numpy.ndarray)
            filename (str)
        """
        driver = gdal.GetDriverByName('GTiff')
        x_size = arr.shape[1]
        y_size = arr.shape[0]
        dataset = driver.Create(filename, x_size, y_size,
                                eType=gdal.GDT_Float32)
        dataset.GetRasterBand(1).WriteArray(arr)