Python gdal.GDT_Float32() Examples

def writeFile(filename,geotransform,geoprojection,data):
	(x,y) = data.shape
	format = "GTiff"
	noDataValue = -9999
	driver = gdal.GetDriverByName(format)
	# you can change the dataformat but be sure to be able to store negative values including -9999
	dst_datatype = gdal.GDT_Float32

	#print(data)

	dst_ds = driver.Create(filename,y,x,1,dst_datatype)
	dst_ds.GetRasterBand(1).WriteArray(data)
	dst_ds.GetRasterBand(1).SetNoDataValue( noDataValue )
	dst_ds.SetGeoTransform(geotransform)
	dst_ds.SetProjection(geoprojection)
	return 1 

def writeFile(filename,geotransform,geoprojection,data):
    (x,y) = data.shape
    format = "GTiff"
    noDataValue = -9999
    driver = gdal.GetDriverByName(format)
    # you can change the dataformat but be sure to be able to store negative values including -9999
    dst_datatype = gdal.GDT_Float32

    #print(data)

    dst_ds = driver.Create(filename,y,x,1,dst_datatype)
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds.GetRasterBand(1).SetNoDataValue( noDataValue )
    dst_ds.SetGeoTransform(geotransform)
    dst_ds.SetProjection(geoprojection)
    return 1 

def calc_ndvi(raster_path, output_path):
    raster = gdal.Open(raster_path)
    out_raster = create_matching_dataset(raster, output_path, datatype=gdal.GDT_Float32)
    array = raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    R = array[2, ...]
    I = array[3, ...]
    out_array[...] = (R-I)/(R+I)

    out_array[...] = np.where(out_array == -2147483648, 0, out_array)

    R = None
    I = None
    array = None
    out_array = None
    raster = None
    out_raster = None 

def save_img(data, geotransform, proj, outPath, noDataValue=np.nan, fieldNames=[]):

    # Start the gdal driver for GeoTIFF
    if outPath == "MEM":
        driver = gdal.GetDriverByName("MEM")
        driverOpt = []

    shape = data.shape
    if len(shape) > 2:
        ds = driver.Create(outPath, shape[1], shape[0], shape[2], gdal.GDT_Float32, driverOpt)
        ds.SetProjection(proj)
        ds.SetGeoTransform(geotransform)
        for i in range(shape[2]):
            ds.GetRasterBand(i+1).WriteArray(data[:, :, i])
            ds.GetRasterBand(i+1).SetNoDataValue(noDataValue)
    else:
        ds = driver.Create(outPath, shape[1], shape[0], 1, gdal.GDT_Float32, driverOpt)
        ds.SetProjection(proj)
        ds.SetGeoTransform(geotransform)
        ds.GetRasterBand(1).WriteArray(data)
        ds.GetRasterBand(1).SetNoDataValue(noDataValue)

    return ds 

def make_geotiff(data, length, width, latn_p, lonw_p, dlat, dlon, outfile, compress_option):

    if data.dtype == np.float32:
        dtype = gdal.GDT_Float32
        nodata = np.nan  ## or 0?
    elif data.dtype == np.uint8:
        dtype = gdal.GDT_Byte
        nodata = None

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(outfile, width, length, 1, dtype, options=compress_option)
    outRaster.SetGeoTransform((lonw_p, dlon, 0, latn_p, 0, dlat))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(data)
    if nodata is not None: outband.SetNoDataValue(nodata)
    outRaster.SetMetadataItem('AREA_OR_POINT', 'Point')
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

    return


#%% Main 

def Save_as_MEM(data='', geo='', projection=''):
    """
    This function save the array as a memory file

    Keyword arguments:
    data -- [array], dataset of the geotiff
    geo -- [minimum lon, pixelsize, rotation, maximum lat, rotation,
            pixelsize], (geospatial dataset)
    projection -- interger, the EPSG code
    """
    # save as a geotiff
    driver = gdal.GetDriverByName("MEM")
    dst_ds = driver.Create('', int(data.shape[1]), int(data.shape[0]), 1,
                           gdal.GDT_Float32, ['COMPRESS=LZW'])
    srse = osr.SpatialReference()
    if projection == '':
        srse.SetWellKnownGeogCS("WGS84")
    else:
        srse.SetWellKnownGeogCS(projection)
    dst_ds.SetProjection(srse.ExportToWkt())
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.SetGeoTransform(geo)
    dst_ds.GetRasterBand(1).WriteArray(data)
    return(dst_ds) 

def calc_raster_terrain_fixed_elevation(crs, elevation, grid_size, raster_path, locator, x_max, x_min, y_max,
                                        y_min):
    # local variables:
    temp_shapefile = locator.get_temporary_file("terrain.shp")
    cols = int((x_max - x_min) / grid_size)
    rows = int((y_max - y_min) / grid_size)
    shapes = Polygon([[x_min, y_min], [x_max, y_min], [x_max, y_max], [x_min, y_max], [x_min, y_min]])
    geodataframe = Gdf(index=[0], crs=crs, geometry=[shapes])
    geodataframe.to_file(temp_shapefile)
    # 1) opening the shapefile
    source_ds = ogr.Open(temp_shapefile)
    source_layer = source_ds.GetLayer()
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Float32)  ##COMMENT 2
    target_ds.SetGeoTransform((x_min, grid_size, 0, y_max, 0, -grid_size))  ##COMMENT 3
    # 5) Adding a spatial reference ##COMMENT 4
    target_dsSRS = osr.SpatialReference()
    target_dsSRS.ImportFromProj4(crs)
    target_ds.SetProjection(target_dsSRS.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(-9999)  ##COMMENT 5
    gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[elevation])  ##COMMENT 6
    target_ds = None  # closing the file 

def write_output_geotiff(md, gt, wkt, data, dest, nodata):
    # pylint: disable=too-many-arguments
    """
    Writes PyRate output data to a GeoTIFF file.

    :param dict md: Dictionary containing PyRate metadata
    :param list gt: GDAL geotransform for the data
    :param list wkt: GDAL projection information for the data
    :param ndarray data: Output data array to save
    :param str dest: Destination file name
    :param float nodata: No data value of data

    :return None, file saved to disk
    """

    driver = gdal.GetDriverByName("GTiff")
    nrows, ncols = data.shape
    ds = driver.Create(dest, ncols, nrows, 1, gdal.GDT_Float32, options=['compress=packbits'])
    # set spatial reference for geotiff
    ds.SetGeoTransform(gt)
    ds.SetProjection(wkt)
    ds.SetMetadataItem(ifc.EPOCH_DATE, str(md[ifc.EPOCH_DATE]))

    # set other metadata
    ds.SetMetadataItem('DATA_TYPE', str(md['DATA_TYPE']))

    # sequence position for time series products

    for k in [ifc.SEQUENCE_POSITION, ifc.PYRATE_REFPIX_X, ifc.PYRATE_REFPIX_Y, ifc.PYRATE_REFPIX_LAT,
              ifc.PYRATE_REFPIX_LON, ifc.PYRATE_MEAN_REF_AREA, ifc.PYRATE_STDDEV_REF_AREA]:
        if k in md:
            ds.SetMetadataItem(k, str(md[k]))

    # write data to geotiff
    band = ds.GetRasterBand(1)
    band.SetNoDataValue(nodata)
    band.WriteArray(data, 0, 0) 

def gdal_dataset(out_fname, columns, rows, driver="GTiff", bands=1,
                 dtype='float32', metadata=None, crs=None,
                 geotransform=None, creation_opts=None):
    """
    Initialises a py-GDAL dataset object for writing image data.
    """
    if dtype == 'float32':
        gdal_dtype = gdal.GDT_Float32
    elif dtype == 'int16':
        gdal_dtype = gdal.GDT_Int16
    else:
        # assume gdal.GDT val is passed to function
        gdal_dtype = dtype

    # create output dataset
    driver = gdal.GetDriverByName(driver)
    outds = driver.Create(out_fname, columns, rows, bands, gdal_dtype, options=creation_opts)

    # geospatial info
    outds.SetGeoTransform(geotransform)
    outds.SetProjection(crs)

    # add metadata
    if metadata is not None:
        for k, v in metadata.items():
            outds.SetMetadataItem(k, str(v))

    return outds 

def Save_as_tiff(name='', data='', geo='', projection=''):
    """
    This function save the array as a geotiff

    Keyword arguments:
    name -- string, directory name
    data -- [array], dataset of the geotiff
    geo -- [minimum lon, pixelsize, rotation, maximum lat, rotation,
            pixelsize], (geospatial dataset)
    projection -- integer, the EPSG code
    """
    # save as a geotiff
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(name, int(data.shape[1]), int(data.shape[0]), 1,
                           gdal.GDT_Float32, ['COMPRESS=LZW'])
    srse = osr.SpatialReference()
    if projection == '':
        srse.SetWellKnownGeogCS("WGS84")

    else:
        try:
            if not srse.SetWellKnownGeogCS(projection) == 6:
                srse.SetWellKnownGeogCS(projection)
            else:
                try:
                    srse.ImportFromEPSG(int(projection))
                except:
                    srse.ImportFromWkt(projection)
        except:
            try:
                srse.ImportFromEPSG(int(projection))
            except:
                srse.ImportFromWkt(projection)

    dst_ds.SetProjection(srse.ExportToWkt())
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.SetGeoTransform(geo)
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds = None
    return() 

def export_raster(self, raster, fname, title, res=None):
        """ Write array to raster file with gdal

        Parameters
        ----------
        raster :    ndarray
                    raster to be exported
        fname :     str
                    file name
        title :     str
                    gdal title of the file
        res :       int/float, optional
                    resolution of the raster in m, if not provided the same as
                    the input CHM
        """
        res = res if res else self.resolution

        fname.parent.mkdir(parents=True, exist_ok=True)

        driver = gdal.GetDriverByName('GTIFF')
        y_pixels, x_pixels = raster.shape
        gdal_file = driver.Create(
            f'{fname}', x_pixels, y_pixels, 1, gdal.GDT_Float32
        )
        gdal_file.SetGeoTransform(
            (self.ul_lon, res, 0., self.ul_lat, 0., -res)
        )
        dataset_srs = gdal.osr.SpatialReference()
        dataset_srs.ImportFromEPSG(self.epsg)
        gdal_file.SetProjection(dataset_srs.ExportToWkt())
        band = gdal_file.GetRasterBand(1)
        band.SetDescription(title)
        band.SetNoDataValue(0.)
        band.WriteArray(raster)
        gdal_file.FlushCache()
        gdal_file = None 

def mk_geotiff_obj(raster, fn, bands=1, gdal_data_type=gdal.GDT_Float32,
                   lat=[46, 45], lon=[-73, -72]):
    """
    Creates a new geotiff file objects using the WGS84 coordinate system, saves
    it to disk, and returns a handle to the python file object and driver

    Parameters
    ------------
    raster : array
        Numpy array of the raster data to be added to the object
    fn : str
        Name of the geotiff file
    bands : int (optional)
        See :py:func:`gdal.GetDriverByName('Gtiff').Create
    gdal_data : gdal.GDT_<type>
        Gdal data type (see gdal.GDT_...)
    lat : list
        northern lat, southern lat
    lon : list
        [western lon, eastern lon]
    """
    NNi, NNj = raster.shape
    driver = gdal.GetDriverByName('GTiff')
    obj = driver.Create(fn, NNj, NNi, bands, gdal_data_type)
    pixel_height = -np.abs(lat[0] - lat[1]) / (NNi - 1.0)
    pixel_width = np.abs(lon[0] - lon[1]) / (NNj - 1.0)
    obj.SetGeoTransform([lon[0], pixel_width, 0, lat[0], 0, pixel_height])
    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS('WGS84')
    obj.SetProjection(srs.ExportToWkt())
    obj.GetRasterBand(1).WriteArray(raster)
    return obj, driver 

def getGDALGDT(dt):
    """
    Need arr.dtype.name in entry.
    Returns gdal_dt from dt (numpy/scipy).

    Parameters
    ----------
    dt : datatype

    Return
    ----------
    gdal_dt : gdal datatype
    """
    if dt == 'bool' or dt == 'uint8':
        gdal_dt = gdal.GDT_Byte
    elif dt == 'int8' or dt == 'int16':
        gdal_dt = gdal.GDT_Int16
    elif dt == 'uint16':
        gdal_dt = gdal.GDT_UInt16
    elif dt == 'int32':
        gdal_dt = gdal.GDT_Int32
    elif dt == 'uint32':
        gdal_dt = gdal.GDT_UInt32
    elif dt == 'int64' or dt == 'uint64' or dt == 'float16' or dt == 'float32':
        gdal_dt = gdal.GDT_Float32
    elif dt == 'float64':
        gdal_dt = gdal.GDT_Float64
    elif dt == 'complex64':
        gdal_dt = gdal.GDT_CFloat64
    else:
        print('Data type non-suported')
        # exit()

    return gdal_dt 

def getDTfromGDAL(gdal_dt):
    """
    Returns datatype (numpy/scipy) from gdal_dt.

    Parameters
    ----------
    gdal_dt : datatype
        data.GetRasterBand(1).DataType

    Return
    ----------
    dt : datatype
    """
    if gdal_dt == gdal.GDT_Byte:
        dt = 'uint8'
    elif gdal_dt == gdal.GDT_Int16:
        dt = 'int16'
    elif gdal_dt == gdal.GDT_UInt16:
        dt = 'uint16'
    elif gdal_dt == gdal.GDT_Int32:
        dt = 'int32'
    elif gdal_dt == gdal.GDT_UInt32:
        dt = 'uint32'
    elif gdal_dt == gdal.GDT_Float32:
        dt = 'float32'
    elif gdal_dt == gdal.GDT_Float64:
        dt = 'float64'
    elif gdal_dt == gdal.GDT_CInt16 or gdal_dt == gdal.GDT_CInt32 or gdal_dt == gdal.GDT_CFloat32 or gdal_dt == gdal.GDT_CFloat64:
        dt = 'complex64'
    else:
        print('Data type unkown')
        # exit()
    return dt 

def reproject_dataset(dataset, output_file_name, wkt_from="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs", epsg_to=32643, pixel_spacing=1000, file_format='GTiff'):
    '''
    :param dataset: Modis subset name use gdal.GetSubdatasets()
    :param output_file_name: file location along with proper extension
    :param wkt_from: Modis wkt (default)
    :param epsg_to: integer default(32643)
    :param pixel_spacing: in metres
    :param file_format: default GTiff
    :return:
    '''
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(epsg_to)
    modis = osr.SpatialReference()
    modis.ImportFromProj4(wkt_from)
    tx = osr.CoordinateTransformation(modis, wgs84)
    g = gdal.Open(dataset)
    geo_t = g.GetGeoTransform()
    print geo_t
    x_size = g.RasterXSize
    y_size = g.RasterYSize
    (ulx, uly, ulz) = tx.TransformPoint(geo_t[0], geo_t[3])
    (lrx, lry, lrz) = tx.TransformPoint(geo_t[0] + (geo_t[1]*x_size), geo_t[3]+ (geo_t[5]*y_size))
    mem_drv = gdal.GetDriverByName(file_format)
    dest = mem_drv.Create(output_file_name, int((lrx-ulx)/pixel_spacing), int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)
    new_geo = ([ulx, pixel_spacing, geo_t[2], uly, geo_t[4], -pixel_spacing])
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(wgs84.ExportToWkt())
    gdal.ReprojectImage(g, dest, modis.ExportToWkt(), wgs84.ExportToWkt(), gdal.GRA_Bilinear)
    print "reprojected" 

def test_band_maths():
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    try:
        os.remove("test_outputs/ndvi.tif")
    except FileNotFoundError:
        pass

    test_file = "test_data/bands.tif"
    test_outputs = "test_outputs/ndvi.tif"
    pyeo.raster_manipulation.apply_band_function(test_file,
                                                 pyeo.raster_manipulation.ndvi_function,
                                                 [2, 3],
                                                 test_outputs,
                                                 out_datatype=gdal.GDT_Float32) 

def mk_test_multifile(testnum, NN, testdir, nx_grid=3, ny_grid=4, nx_overlap=16,
                      ny_overlap=32, lat=[46, 45], lon=[-73, -72]):
    """
    Written to make test case for multi-file edge resolution.
    """
    path = os.path.split(make_file_names(testnum, NN, os.path.join(
                                               testdir, 'chunks'))['elev'])[0]
    try:
        os.makedirs(path)
    except:
        pass

    def _get_chunk_edges(NN, chunk_size, chunk_overlap):
        left_edge = np.arange(0, NN - chunk_overlap, chunk_size)
        left_edge[1:] -= chunk_overlap
        right_edge = np.arange(0, NN - chunk_overlap, chunk_size)
        right_edge[:-1] = right_edge[1:] + chunk_overlap
        right_edge[-1] = NN
        right_edge = np.minimum(right_edge, NN)
        return left_edge, right_edge

    elev_data, ang_data, fel_data = get_test_data(testnum, NN)
    try:
        raster = fel_data.raster_data
    except:
        raster = elev_data.raster_data

    ni, nj = raster.shape

    top_edge, bottom_edge = _get_chunk_edges(ni, ni // ny_grid, ny_overlap)
    left_edge, right_edge = _get_chunk_edges(nj, nj // nx_grid, nx_overlap)

#    gc = elev_data.grid_coordinates
#    lat = gc.y_axis
#    lon = gc.x_axis
    lat = np.linspace(lat[0], lat[1], ni)
    lon = np.linspace(lon[0], lon[1], nj)
    count = 0
    for te, be in zip(top_edge, bottom_edge):
        for le, re in zip(left_edge, right_edge):
            count += 1
            fn = os.path.join(path,
                              get_fn_from_coords((lat[be-1], lon[le], lat[te],
                                                  lon[re-1]), 'elev'))
            print count, ": [%d:%d, %d:%d]" % (te, be, le, re), \
                '(lat, lon) = (%g to %g, %g to %g)' % (lat[te], lat[be-1],
                                                       lon[le], lon[re-1]), \
                'min,max = (%g to %g, %g to %g)' % (lat.min(), lat.max(),
                                                    lon.min(), lon.max())
            mk_geotiff_obj(raster[te:be, le:re], fn,
                           bands=1, gdal_data_type=gdal.GDT_Float32,
                           lat=[lat[te], lat[be-1]], lon=[lon[le], lon[re-1]])


# %% MAKE ALL THE TESTS 

def Vector_to_Raster(Dir, shapefile_name, reference_raster_data_name):
    """
    This function creates a raster of a shp file

    Keyword arguments:
    Dir --
        str: path to the basin folder
    shapefile_name -- 'C:/....../.shp'
        str: Path from the shape file
    reference_raster_data_name -- 'C:/....../.tif'
        str: Path to an example tiff file (all arrays will be reprojected to this example)
    """

    from osgeo import gdal, ogr

    geo, proj, size_X, size_Y=Open_array_info(reference_raster_data_name)

    x_min = geo[0]
    x_max = geo[0] + size_X * geo[1]
    y_min = geo[3] + size_Y * geo[5]
    y_max = geo[3]
    pixel_size = geo[1]

    # Filename of the raster Tiff that will be created
    Dir_Basin_Shape = os.path.join(Dir,'Basin')
    if not os.path.exists(Dir_Basin_Shape):
        os.mkdir(Dir_Basin_Shape)

    Basename = os.path.basename(shapefile_name)
    Dir_Raster_end = os.path.join(Dir_Basin_Shape, os.path.splitext(Basename)[0]+'.tif')

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

    # Create the destination data source
    x_res = int(round((x_max - x_min) / pixel_size))
    y_res = int(round((y_max - y_min) / pixel_size))

    # Create tiff file
    target_ds = gdal.GetDriverByName('GTiff').Create(Dir_Raster_end, x_res, y_res, 1, gdal.GDT_Float32, ['COMPRESS=LZW'])
    target_ds.SetGeoTransform(geo)
    srse = osr.SpatialReference()
    srse.SetWellKnownGeogCS(proj)
    target_ds.SetProjection(srse.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    target_ds.GetRasterBand(1).SetNoDataValue(-9999)
    band.Fill(-9999)

    # Rasterize the shape and save it as band in tiff file
    gdal.RasterizeLayer(target_ds, [1], source_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
    target_ds = None

    # Open array
    Raster_Basin = Open_tiff_array(Dir_Raster_end)

    return(Raster_Basin) 

def Transform(SourceLANDSAF, OutPath, Type, Dates = ['2000-01-01','2013-12-31']):
    """
    This function creates short wave maps based on the SIS and SID 
    This function converts packed nc files to gtiff file of comparable file size
  		
    Keyword arguments:
    SourceLANDSAF -- 'C:/'  path to the LANDSAF source data (The directory includes SIS and SID)
    Dir -- 'C:/' path to the WA map
    latlim -- [ymin, ymax] (values must be between -60 and 60)
    lonlim -- [xmin, xmax] (values must be between -180 and 180)
    Dates -- ['yyyy-mm-dd','yyyy-mm-dd']
    """
			
    path = os.path.join(SourceLANDSAF,Type)

    os.chdir(path)
    Dates = pd.date_range(Dates[0],Dates[1],freq='D')

    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS("WGS84")  
    projection = srs.ExportToWkt() 
    driver = gdal.GetDriverByName("GTiff")
    
    
    for Date in Dates:
        if Type == 'SIS': 
            ZipFile = glob.glob('SISdm%s*.nc.gz' % Date.strftime('%Y%m%d'))[0]		
            File = os.path.splitext(ZipFile)[0] 
        elif Type == 'SID':
            ZipFile = glob.glob('*dm%s*.nc.gz' % Date.strftime('%Y%m%d'))[0]		
            File = os.path.splitext(ZipFile)[0]
        
        # find path to the executable 	 
        fullCmd = ''.join("7z x %s -o%s -aoa"  %(os.path.join(path,ZipFile),path))   
        process = subprocess.Popen(fullCmd)
        process.wait()
	
        NC = netCDF4.Dataset(File,'r+',format='NETCDF4')
        Data = NC[Type][0,:,:]
        lon = NC.variables['lon'][:][0]	- 0.025	
        lat = NC.variables['lat'][:][-1] + 0.025					
        geotransform = [lon,0.05,0,lat,0,-0.05]								
								
        dst_ds = driver.Create(OutPath, int(np.size(Data,1)), int(np.size(Data,0)),  1, gdal.GDT_Float32,  ['COMPRESS=DEFLATE'])
        # set the reference info 
        dst_ds.SetProjection(projection)
        dst_ds.SetGeoTransform(geotransform)
        dst_ds.GetRasterBand(1).SetNoDataValue(-1)
        dst_ds.GetRasterBand(1).WriteArray(np.flipud(Data))
        NC.close()        
        del dst_ds, NC, Data
        
        os.remove(File) 

def open_data(filename):
    '''
    The function open and load the image given its name.
    The type of the data is checked from the file and the scipy array is initialized accordingly.
    Input:
        filename: the name of the file
    Output:
        im: the data cube
        GeoTransform: the geotransform information
        Projection: the projection information
    '''
    data = gdal.Open(filename, gdal.GA_ReadOnly)
    if data is None:
        print('Impossible to open ' + filename)
        # exit()
    nc = data.RasterXSize
    nl = data.RasterYSize
    d = data.RasterCount

    # Get the type of the data
    gdal_dt = data.GetRasterBand(1).DataType
    if gdal_dt == gdal.GDT_Byte:
        dt = 'uint8'
    elif gdal_dt == gdal.GDT_Int16:
        dt = 'int16'
    elif gdal_dt == gdal.GDT_UInt16:
        dt = 'uint16'
    elif gdal_dt == gdal.GDT_Int32:
        dt = 'int32'
    elif gdal_dt == gdal.GDT_UInt32:
        dt = 'uint32'

    elif gdal_dt == gdal.GDT_Float32:
        dt = 'float32'
    elif gdal_dt == gdal.GDT_Float64:
        dt = 'float64'
    elif gdal_dt == gdal.GDT_CInt16 or gdal_dt == gdal.GDT_CInt32 or gdal_dt == gdal.GDT_CFloat32 or gdal_dt == gdal.GDT_CFloat64:
        dt = 'complex64'
    else:
        print('Data type unkown')
        # exit()

    # Initialize the array
    if d == 1:
        im = np.empty((nl, nc), dtype=dt)
    else:
        im = np.empty((nl, nc, d), dtype=dt)

    if d == 1:
        im[:, :] = data.GetRasterBand(1).ReadAsArray()
    else:
        for i in range(d):
            im[:, :, i] = data.GetRasterBand(i + 1).ReadAsArray()

    GeoTransform = data.GetGeoTransform()
    Projection = data.GetProjection()
    data = None
    return im, GeoTransform, Projection