Python gdal.GDT_Int32() Examples

def save_raster_simple(array, path, dst_filename):

    """ Save an array base on an existing raster """

    example = gdal.Open(path)
    x_pixels = array.shape[1]  # number of pixels in x
    y_pixels = array.shape[0]  # number of pixels in y
    bands = 1
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(dst_filename,x_pixels, y_pixels, bands, 
                            gdal.GDT_Int32)

    geotrans=example.GetGeoTransform()  #get GeoTranform from existed 'data0'
    proj=example.GetProjection() #you can get from a exsited tif or import 
    dataset.SetGeoTransform(geotrans)
    dataset.SetProjection(proj)

    dataset.GetRasterBand(1).WriteArray(array[:,:])

    dataset.FlushCache() 

def save_raster(array, path, dst_filename):
    """ Save the final multiband array based on an existing raster """

    example = gdal.Open(path)
    x_pixels = array.shape[2]  # number of pixels in x
    y_pixels = array.shape[1]  # number of pixels in y
    bands = array.shape[0]
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(dst_filename,x_pixels, 
                            y_pixels, bands ,gdal.GDT_Int32)

    geotrans=example.GetGeoTransform()  #get GeoTranform from existed 'data0'
    proj=example.GetProjection() #you can get from a exsited tif or import 
    dataset.SetGeoTransform(geotrans)
    dataset.SetProjection(proj)

    for b in range(bands):
        dataset.GetRasterBand(b+1).WriteArray(array[b,:,:])

    dataset.FlushCache() 

def save_raster_memory(array, path):

    """ Save a raster into memory """

    example = gdal.Open(path)
    x_pixels = array.shape[1]  # number of pixels in x
    y_pixels = array.shape[0]  # number of pixels in y
    driver = gdal.GetDriverByName('MEM')
    dataset = driver.Create('',x_pixels, y_pixels, 1,gdal.GDT_Int32)
    dataset.GetRasterBand(1).WriteArray(array[:,:])

    # follow code is adding GeoTranform and Projection
    geotrans=example.GetGeoTransform()  #get GeoTranform from existed 'data0'
    proj=example.GetProjection() #you can get from a exsited tif or import 
    dataset.SetGeoTransform(geotrans)
    dataset.SetProjection(proj)

    return dataset 

def apply_image_function(in_paths, out_path, function, out_datatype = gdal.GDT_Int32):
    """Applies a pixel-wise function across every image. Assumes each image is exactly contiguous and, for now,
    single-banded. function() should take a list of values and return a single value."""
    rasters = [gdal.Open(in_path) for in_path in in_paths]
    raster_arrays = [raster.GetVirtualMemArray() for raster in rasters]
    in_array = np.stack(raster_arrays, axis=0)

    out_raster = create_matching_dataset(rasters[0], out_path=out_path, datatype=out_datatype)
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    out_array[...] = np.apply_along_axis(function, 0, in_array)

    # Deallocating. Not taking any chances here.
    out_array = None
    out_raster = None
    in_array = None
    for raster_array, raster in zip(raster_arrays, rasters):
        raster_array = None
        raster = None 

def apply_band_function(in_path, function, bands, out_path, out_datatype = gdal.GDT_Int32):
    """Applys an arbitrary band mathemtics function to an image at in_path and saves the result at out_map.
    Function should be a function ofblect of the form f(band_input_A, band_input_B, ...)"""
    raster = gdal.Open(in_path)
    out_raster = create_matching_dataset(raster, out_path=out_path, datatype=out_datatype)
    array = raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    band_views = [array[band, ...] for band in bands]
    out_array[...] = function(*band_views)
    out_array = None
    for view in band_views:
        view = None
    raster = None
    out_raster = None 

def save_tif(self, mask, file_path):
        '''Save the given mask as a .tif file.

        Parameters
            mask 		-	A mask generated with masker.
            file_path	-	Path of .tif file.
        '''
        import gdal

        driver = gdal.GetDriverByName('GTiff')

        x_pixels = mask.shape[1]
        y_pixels = mask.shape[0]

        dataset = driver.Create(file_path, x_pixels, y_pixels, 1, gdal.GDT_Int32)

        if self.file_path is not None:
            extension = os.path.splitext(self.file_path)[1].lower()
            if extension == '.hdf':
                hdfdataset = gdal.Open(self.file_path)
                subdataset = hdfdataset.GetSubDatasets()[0][0]
                bandfile = gdal.Open(subdataset)
            else:
                bandfile = gdal.Open(self.file_path)

            dataset.SetGeoTransform(bandfile.GetGeoTransform())
            dataset.SetProjection(bandfile.GetProjectionRef())

        dataset.GetRasterBand(1).WriteArray(mask)
        dataset.FlushCache() 

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 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 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 mosaic_images(raster_paths, out_raster_file, format="GTiff", datatype=gdal.GDT_Int32, nodata = 0):
    """
    Mosaics multiple images with the same number of layers into one single image. Overwrites
    overlapping pixels with the value furthest down raster_paths. Takes projection from the first
    raster.

    Parameters
    ----------
    raster_paths
        A list of paths of raster to be mosaiced
    out_raster_file
        The path to the output file
    format
        The image format of the output raster.
    datatype
        The datatype of the output raster
    nodata
        The input nodata value; any pixels in raster_paths with this value will be ignored.

    """

    # This, again, is very similar to stack_rasters
    log = logging.getLogger(__name__)
    log.info("Beginning mosaic")
    rasters = [gdal.Open(raster_path) for raster_path in raster_paths]
    projection = rasters[0].GetProjection()
    in_gt = rasters[0].GetGeoTransform()
    x_res = in_gt[1]
    y_res = in_gt[5] * -1  # Y resolution in agt is -ve for Maths reasons
    combined_polygon = align_bounds_to_whole_number(get_combined_polygon(rasters, geometry_mode='union'))
    layers = rasters[0].RasterCount
    out_raster = create_new_image_from_polygon(combined_polygon, out_raster_file, x_res, y_res, layers,
                                               projection, format, datatype)
    log.info("New empty image created at {}".format(out_raster_file))
    out_raster_array = out_raster.GetVirtualMemArray(eAccess=gdal.GF_Write)
    for i, raster in enumerate(rasters):
        log.info("Now mosaicking raster no. {}".format(i))
        in_raster_array = raster.GetVirtualMemArray()
        if len(in_raster_array.shape) == 2:
            in_raster_array = np.expand_dims(in_raster_array, 0)
        in_bounds = get_raster_bounds(raster)
        out_x_min, out_x_max, out_y_min, out_y_max = pixel_bounds_from_polygon(out_raster, in_bounds)
        out_raster_view = out_raster_array[:, out_y_min: out_y_max, out_x_min: out_x_max]
        np.copyto(out_raster_view, in_raster_array, where=in_raster_array != nodata)
        in_raster_array = None
        out_raster_view = None
    log.info("Raster mosaicking done")
    out_raster_array = None 

def create_new_image_from_polygon(polygon, out_path, x_res, y_res, bands,
                           projection, format="GTiff", datatype = gdal.GDT_Int32, nodata = -9999):
    """
    Returns an empty image that covers the extent of the imput polygon.

    Parameters
    ----------
    polygon
        An OGR.Geometry object of a single polygon
    out_path
        The path to save the new image to
    x_res
        Pixel width in the new image
    y_res
        Pixel height in the new image
    bands
        Number of bands in the new image.
    projection
        The projection, in wkt, of the output image.
    format
        The gdal raster format of the output image
    datatype
        The gdal datatype of the output image
    nodata
        The nodata value of the output image

    Returns
    -------
    A gdal.Image object

    """
    # TODO: Implement nodata
    bounds_x_min, bounds_x_max, bounds_y_min, bounds_y_max = polygon.GetEnvelope()
    if bounds_x_min >= bounds_x_max:
        bounds_x_min, bounds_x_max = bounds_x_max, bounds_x_min
    if bounds_y_min >= bounds_y_max:
        bounds_y_min, bounds_y_max = bounds_y_max, bounds_y_min
    final_width_pixels = int(np.abs(bounds_x_max - bounds_x_min) / x_res)
    final_height_pixels = int(np.abs(bounds_y_max - bounds_y_min) / y_res)
    driver = gdal.GetDriverByName(format)
    out_raster = driver.Create(
        out_path, xsize=final_width_pixels, ysize=final_height_pixels,
        bands=bands, eType=datatype
    )
    out_raster.SetGeoTransform([
        bounds_x_min, x_res, 0,
        bounds_y_max, 0, y_res * -1
    ])
    out_raster.SetProjection(projection)
    return out_raster 

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