Python osgeo.gdal.GDT_Float32() Examples

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 _create_geotiff_file(outfn, driver, shape, metadata, num_bands):
    dst = driver.Create(
        outfn,
        shape[1],
        shape[0],
        num_bands,
        gdal.GDT_Float32,
        ["COMPRESS=DEFLATE", "PREDICTOR=3"],
    )

    sx = (metadata["x2"] - metadata["x1"]) / shape[1]
    sy = (metadata["y2"] - metadata["y1"]) / shape[0]
    dst.SetGeoTransform([metadata["x1"], sx, 0.0, metadata["y2"], 0.0, -sy])

    sr = osr.SpatialReference()
    sr.ImportFromProj4(metadata["projection"])
    dst.SetProjection(sr.ExportToWkt())

    return dst 

def save_geotiff(self,filename, dtype=gdal.GDT_Float32, **kwargs):
        """
        Save a GeoTIFF of the raster currently loaded.

        Georeferenced and subset according to the current raster.

        :params filename: the path and filename to save the file to.
        :type filename: str
        :params dtype: GDAL datatype, defaults to Float32.
        :type dtype: int

        :returns: 1 on success.
        :rtype: int

        """

        if self.r is None:
            raise ValueError('There are no image data loaded. No file can be created.')

        simple_write_geotiff(filename, self.r, self.trans, 
            wkt=self.srs.ExportToWkt(), dtype=dtype, **kwargs) 

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 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 test_gdal_create_dataset(self):
        testfunc = io.gdal.gdal_create_dataset
        tmp = tempfile.NamedTemporaryFile(mode="w+b").name
        with pytest.raises(TypeError):
            testfunc("AIG", tmp)
        from osgeo import gdal

        with pytest.raises(TypeError):
            testfunc(
                "AAIGrid", tmp, cols=10, rows=10, bands=1, gdal_type=gdal.GDT_Float32
            )
        testfunc("GTiff", tmp, cols=10, rows=10, bands=1, gdal_type=gdal.GDT_Float32)
        testfunc(
            "GTiff",
            tmp,
            cols=10,
            rows=10,
            bands=1,
            gdal_type=gdal.GDT_Float32,
            remove=True,
        ) 

def getNumpyType(self, pixelType = gdal.GDT_Byte):
        """
        Translates the gdal raster type to numpy type
        pixelType: gdal raster type
        """
        if pixelType == gdal.GDT_Byte:
            return numpy.uint8
        elif pixelType == gdal.GDT_UInt16:
            return numpy.uint16
        elif pixelType == gdal.GDT_Int16:
            return numpy.int16
        elif pixelType == gdal.GDT_UInt32:
            return numpy.uint32
        elif pixelType == gdal.GDT_Int32:
            return numpy.int32
        elif pixelType == gdal.GDT_Float32:
            return numpy.float32
        elif pixelType == gdal.GDT_Float64:
            return numpy.float64 

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 getNumpyType(self, pixelType = gdal.GDT_Byte):
        '''
        Translates the gdal raster type to numpy type
        pixelType: gdal raster type
        '''
        if pixelType == gdal.GDT_Byte:
            return numpy.uint8
        elif pixelType == gdal.GDT_UInt16:
            return numpy.uint16
        elif pixelType == gdal.GDT_Int16:
            return numpy.int16
        elif pixelType == gdal.GDT_UInt32:
            return numpy.uint32
        elif pixelType == gdal.GDT_Int32:
            return numpy.int32
        elif pixelType == gdal.GDT_Float32:
            return numpy.float32
        elif pixelType == gdal.GDT_Float64:
            return numpy.float64 

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 _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 writeGeotif(lat, lon, res, data, filename=None):
    """Writes Geotif in temporary directory so it can be imported into the PostGIS database."""
    if isinstance(data, np.ma.masked_array):
        nodata = np.double(data.fill_value)
        data = data.data
    else:
        nodata = -9999.
    if len(data.shape) < 2:
        nrows = int((max(lat) - min(lat)) / res) + 1
        ncols = int((max(lon) - min(lon)) / res) + 1
        out = np.zeros((nrows, ncols)) + nodata
        for c in range(len(lat)):
            i = int((max(lat) - lat[c]) / res)
            j = int((lon[c] - min(lon)) / res)
            out[i, j] = data[c]
    else:
        nrows, ncols = data.shape
        out = data
    if filename is None:
        f = tempfile.NamedTemporaryFile(suffix=".tif", delete=False)
        filename = f.name
        f.close()
    driver = gdal.GetDriverByName("GTiff")
    ods = driver.Create(filename, ncols, nrows, 1, gdal.GDT_Float32)
    ods.SetGeoTransform([min(lon) - res / 2.0, res, 0,
                         max(lat) + res / 2.0, 0, -res])
    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS("WGS84")
    ods.SetProjection(srs.ExportToWkt())
    ods.GetRasterBand(1).WriteArray(out)
    ods.GetRasterBand(1).SetNoDataValue(nodata)
    ods = None
    return filename 

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

def dump_raster(rast, rast_path, driver='GTiff', gdt=None, nodata=None):
    '''
    Creates a raster file from a given gdal.Dataset instance. Arguments:
        rast        A gdal.Dataset; does NOT accept NumPy array
        rast_path   The path of the output raster file
        driver      The name of the GDAL driver to use (determines file type)
        gdt         The GDAL data type to use, e.g., see gdal.GDT_Float32
        nodata      The NoData value; defaults to -9999.
    '''
    if gdt is None:
        gdt = rast.GetRasterBand(1).DataType
    driver = gdal.GetDriverByName(driver)
    sink = driver.Create(
        rast_path, rast.RasterXSize, rast.RasterYSize, rast.RasterCount, int(gdt))
    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 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])
    raster_field = ogr.FieldDefn('id', type_mapping[gdal.GDT_Int32])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# extract sinks from dem 

def array2raster(rasterfn,newRasterfn,array,driver_name = "ENVI", noDataValue = -9999):
    """Takes an array and writes to a GDAL compatible raster. It needs another raster to map the dimensions.

    Args:
        FileName (str): The filename (with path and extension) of a raster that has the same dimensions as the raster to be written.
        newRasterfn (str): The filename (with path and extension) of the new raster.
        array (np.array): The array to be written
        driver_name (str): The type of raster to write. Default is ENVI since that is the LSDTOpoTools format
        noDataValue (float): The no data value

    Return:
        np.array: A numpy array with the data from the raster.

    Author: SMM
    """

    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName(driver_name)
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outRaster.GetRasterBand(1).SetNoDataValue( noDataValue )
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
#============================================================================== 

def _str_of_gdt(gdt):
    return {
        gdal.GDT_Byte: 'GDT_Byte',
        gdal.GDT_Int16: 'GDT_Int16',
        gdal.GDT_Int32: 'GDT_Int32',
        gdal.GDT_UInt16: 'GDT_UInt16',
        gdal.GDT_UInt32: 'GDT_UInt32',
        gdal.GDT_Float32: 'GDT_Float32',
        gdal.GDT_Float64: 'GDT_Float64',
        gdal.GDT_CFloat32: 'GDT_CFloat32',
        gdal.GDT_CFloat64: 'GDT_CFloat64',
        gdal.GDT_CInt16: 'GDT_CInt16',
        gdal.GDT_CInt32: 'GDT_CInt32',
    }[gdt] 

def make_tif2(path, reso=(1., -1.), rsize=(10, 10), tl=(0., 10.),
              proj=SRS[0]['wkt'], channel_count=1, dtype=gdal.GDT_Float32,
              nodata=-32000, nodata_border_size=(0, 0, 0, 0)):
    """Create a tiff files"""
    reso = np.asarray(reso)
    fp = buzz.Footprint(tl=tl, rsize=rsize, size=np.abs(reso * rsize))
    x, y = fp.meshgrid_raster
    a = x + y
    if nodata_border_size != 0:
        l, r, t, b = nodata_border_size
        if t != 0:
            a[None:t, None:None] = nodata
        if b != 0:
            a[-b:None, None:None] = nodata
        if l != 0:
            a[None:None, None:l] = nodata
        if r != 0:
            a[None:None, -r:None] = nodata

    LOGGER.info('TIFF ARRAY:%s\n', a)
    gdal.UseExceptions()
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(path, int(rsize[0]), int(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(nodata)
    dataset.FlushCache() 

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 gdal_create_dataset(
    drv, name, cols=0, rows=0, bands=0, gdal_type=gdal.GDT_Unknown, remove=False
):
    """Creates GDAL.DataSet object.

    Parameters
    ----------
    drv : string
        GDAL driver string
    name : string
        path to filename
    cols : int
        # of columns
    rows : int
        # of rows
    bands : int
        # of raster bands
    gdal_type : raster data type
        eg. gdal.GDT_Float32
    remove : bool
        if True, existing gdal.Dataset will be
        removed before creation

    Returns
    -------
    out : gdal.Dataset
        object

    """
    driver = gdal.GetDriverByName(drv)
    metadata = driver.GetMetadata()

    if not metadata.get("DCAP_CREATE", False):
        raise TypeError(
            "WRADLIB: Driver {} doesn't support " "Create() method.".format(drv)
        )

    if remove:
        if os.path.exists(name):
            driver.Delete(name)
    ds = driver.Create(name, cols, rows, bands, gdal_type)

    return ds 

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 RasterDifference(RasterFile1, RasterFile2, raster_band=1, OutFileName="Test.outfile", OutFileType="ENVI"):
    """
    Takes two rasters of same size and subtracts second from first,
    e.g. Raster1 - Raster2 = raster_of_difference
    then writes it out to file
    """

    Raster1 = gdal.Open(RasterFile1)
    Raster2 = gdal.Open(RasterFile2)

    print("RASTER 1: ")
    print(Raster1.GetGeoTransform())
    print(Raster1.RasterCount)
    print(Raster1.GetRasterBand(1).XSize)
    print(Raster1.GetRasterBand(1).YSize)
    print(Raster1.GetRasterBand(1).DataType)

    print("RASTER 2: ")
    print(Raster2.GetGeoTransform())
    print(Raster2.RasterCount)
    print(Raster2.GetRasterBand(1).XSize)
    print(Raster2.GetRasterBand(1).YSize)
    print(Raster2.GetRasterBand(1).DataType)

    raster_array1 = np.array(Raster1.GetRasterBand(raster_band).ReadAsArray())
    raster_array2 = np.array(Raster2.GetRasterBand(raster_band).ReadAsArray())

    assert(raster_array1.shape == raster_array2.shape )
    print("Shapes: ", raster_array1.shape, raster_array2.shape)

    difference_raster_array = raster_array1 - raster_array2

#    import matplotlib.pyplot as plt
#
#    plt.imshow(difference_raster_array)
#

    driver = gdal.GetDriverByName(OutFileType)

    dsOut = driver.Create(OutFileName,
                          Raster1.GetRasterBand(1).XSize,
                          Raster1.GetRasterBand(1).YSize,
                          1,
                          gdal.GDT_Float32)
                          #Raster1.GetRasterBand(raster_band).DataType)


    gdal_array.CopyDatasetInfo(Raster1,dsOut)
    bandOut = dsOut.GetRasterBand(1)
    gdal_array.BandWriteArray(bandOut, difference_raster_array)

#============================================================================== 

def _export_geotiff(F, exporter):
    def init_band(band):
        band.SetScale(1.0)
        band.SetOffset(0.0)
        band.SetUnitType(exporter["metadata"]["unit"])

    if exporter["incremental"] is None:
        for i in range(exporter["num_timesteps"]):
            if exporter["num_ens_members"] == 1:
                band = exporter["dst"][i].GetRasterBand(1)
                init_band(band)
                band.WriteArray(F[i, :, :])
            else:
                for j in range(exporter["num_ens_members"]):
                    band = exporter["dst"][i].GetRasterBand(j + 1)
                    init_band(band)
                    band.WriteArray(F[j, i, :, :])
    elif exporter["incremental"] == "timestep":
        i = exporter["num_files_written"]

        outfn = _get_geotiff_filename(
            exporter["outfnprefix"],
            exporter["startdate"],
            exporter["num_timesteps"],
            exporter["timestep"],
            i,
        )
        dst = _create_geotiff_file(
            outfn,
            exporter["driver"],
            exporter["shape"],
            exporter["metadata"],
            exporter["num_ens_members"],
        )

        for j in range(exporter["num_ens_members"]):
            band = dst.GetRasterBand(j + 1)
            init_band(band)
            if exporter["num_ens_members"] > 1:
                band.WriteArray(F[j, :, :])
            else:
                band.WriteArray(F)

        exporter["num_files_written"] += 1
    elif exporter["incremental"] == "member":
        for i in range(exporter["num_timesteps"]):
            # NOTE: This does not work because the GeoTIFF driver does not
            # support adding bands. An alternative solution needs to be
            # implemented.
            exporter["dst"][i].AddBand(gdal.GDT_Float32)
            band = exporter["dst"][i].GetRasterBand(exporter["dst"][i].RasterCount)
            init_band(band)
            band.WriteArray(F[i, :, :])