Python gdal.GA_ReadOnly() Examples

def do_proximity(config, image, srcarray, label, _class):

    """ Get proximity of each pixel to to certain value """

    #First create a band in memory that's that's just 1s and 0s
    src_ds = gdal.Open( image, gdal.GA_ReadOnly )
    srcband = src_ds.GetRasterBand(1)
    mem_rast = save_raster_memory(srcarray, image)
    mem_band = mem_rast.GetRasterBand(1)

    #Now the code behind gdal_sieve.py
    maskband = None
    drv = gdal.GetDriverByName('GTiff')

    #Need to output sieved file
    dst_path = config['postprocessing']['deg_class']['prox_dir']
    
    dst_filename = dst_path + '/' + image.split('/')[-1].split('.')[0] + '_' + _class + '.tif'
    dst_ds = drv.Create( dst_filename,src_ds.RasterXSize, src_ds.RasterYSize,1,
                         srcband.DataType )

    wkt = src_ds.GetProjection()
    if wkt != '':
        dst_ds.SetProjection( wkt )
    dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )

    dstband = dst_ds.GetRasterBand(1)

    # Parameters
    prog_func = None

    options = []
    options.append( 'VALUES=' + str(label) )

    result = gdal.ComputeProximity(mem_band, dstband, options)#, 
                              #callback = prog_func )

    proximity = dstband.ReadAsArray()

    return proximity 

def open_inumpyut(self,vrt):
        if vrt:
            self.in_ds = gdal.Open(vrt, gdal.GA_ReadOnly)
        else:
            raise Exception("No vrt file was specified")
            
        if self.options.verbose:
            print("Inumpyut file:", "( %sP x %sL - %s bands)" % (self.in_ds.RasterXSize, self.in_ds.RasterYSize, self.in_ds.RasterCount))

        if not self.in_ds:
            # Note: GDAL prints the ERROR message too
            self.error("It is not possible to open the inumpyut file '%s'." % vrt )
            
        if self.in_ds.RasterCount == 0:
            self.error( "Inumpyut file '%s' has no raster band" % vrt )
        
        self.out_ds = self.in_ds
        
        # Get alpha band (either directly or from NODATA value)
        self.alphaband = self.out_ds.GetRasterBand(1).GetMaskBand()
        if (self.alphaband.GetMaskFlags() & gdal.GMF_ALPHA) or self.out_ds.RasterCount==4 or self.out_ds.RasterCount==2:
            self.dataBandsCount = self.out_ds.RasterCount - 1
        else:
            self.dataBandsCount = self.out_ds.RasterCount

    # ------------------------------------------------------------------------- 

def sieve_strata(config, image):

    """ First create a band in memory that's that's just 1s and 0s """

    src_ds = gdal.Open( image, gdal.GA_ReadOnly )
    srcband = src_ds.GetRasterBand(1)
    srcarray = srcband.ReadAsArray()

    mem_rast = save_raster_memory(srcarray, image)
    mem_band = mem_rast.GetRasterBand(1)

    #Now the code behind gdal_sieve.py
    maskband = None
    drv = gdal.GetDriverByName('GTiff')

    #Need to output sieved file
    dst_path = config['postprocessing']['sieve']['sieve_file']
    dst_filename = dst_path + '/' + image.split('/')[-1].split('.')[0] + '_sieve_mask.tif'

    dst_ds = drv.Create( dst_filename,src_ds.RasterXSize, src_ds.RasterYSize,1,
                         srcband.DataType )
    wkt = src_ds.GetProjection()
    if wkt != '':
        dst_ds.SetProjection( wkt )
    dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )

    dstband = dst_ds.GetRasterBand(1)

    # Parameters
    prog_func = None

    threshold = config['postprocessing']['sieve']['threshold_strata']
    connectedness = config['postprocessing']['sieve']['connectedness']
    if connectedness not in [8, 4]:
	print "connectness only supports value of 4 or 8"
	sys.exit()
    result = gdal.SieveFilter(mem_band, maskband, dstband,
                              threshold, connectedness,
                              callback = prog_func )

    sieved = dstband.ReadAsArray()

    dst_full = config['postprocessing']['sieve']['sieved_output']

    if dst_full:
	dst_full_name = dst_full + '/' + image.split('/')[-1].split('.')[0] + '_sieved.tif'
        save_raster_simple(sieved, image, dst_full_name)

    return sieved 

def sieve(config, image):

    """ First create a band in memory that's that's just 1s and 0s """

    src_ds = gdal.Open( image, gdal.GA_ReadOnly )
    srcband = src_ds.GetRasterBand(1)
    srcarray = srcband.ReadAsArray()
    srcarray[srcarray > 0] = 1

    mem_rast = save_raster_memory(srcarray, image)
    mem_band = mem_rast.GetRasterBand(1)

    #Now the code behind gdal_sieve.py
    maskband = None
    drv = gdal.GetDriverByName('GTiff')

    #Need to output sieved file
    dst_path = config['postprocessing']['sieve']['sieve_file']
    dst_filename = dst_path + '/' + image.split('/')[-1].split('.')[0] + '_sieve_mask.tif'

    dst_ds = drv.Create( dst_filename,src_ds.RasterXSize, src_ds.RasterYSize,1,
                         srcband.DataType )
    wkt = src_ds.GetProjection()
    if wkt != '':
        dst_ds.SetProjection( wkt )
    dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )

    dstband = dst_ds.GetRasterBand(1)

    # Parameters
    prog_func = None

    threshold = config['postprocessing']['sieve']['threshold']
    connectedness = config['postprocessing']['sieve']['connectedness']
    if connectedness not in [8, 4]:
	print "connectness only supports value of 4 or 8"
	sys.exit()
    result = gdal.SieveFilter(mem_band, maskband, dstband,
                              threshold, connectedness,
                              callback = prog_func )

    sieved = dstband.ReadAsArray()
    sieved[sieved < 0] = 0

    src_new = gdal.Open(image)
    out_img = src_new.ReadAsArray().astype(np.float)

    out_img[np.isnan(out_img)] = 0 

    for b in range(out_img.shape[0]):
        out_img[b,:,:][sieved == 0] = 0

    dst_full = config['postprocessing']['sieve']['sieved_output']

    if dst_full:
	dst_full_name = dst_full + '/' + image.split('/')[-1].split('.')[0] + '_sieved.tif'
        save_raster(out_img, image, dst_full_name)

    return out_img 

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 

def get_metadata_from_loc(loc_file: str, lut_params: LUTConfig, trim_lines: int = 5, nodata_value: float = -9999) -> \
        (float, float, float, np.array):
    """ Get metadata needed for complete runs from the location file (bands long, lat, elev).

    Args:
        loc_file: file name to pull data from
        lut_params: parameters to use to define lut grid
        trim_lines: number of lines to ignore at beginning and end of file (good if lines contain values that are
                    erroneous but not nodata
        nodata_value: value to ignore from location file

    :Returns:
        tuple containing:
            mean_latitude - mean latitude of good values from the location file
            mean_longitude - mean latitude of good values from the location file
            mean_elevation_km - mean ground estimate of good values from the location file
            elevation_lut_grid - the elevation look up table, based on globals and values from location file
    """

    loc_dataset = gdal.Open(loc_file, gdal.GA_ReadOnly)

    loc_data = np.zeros((loc_dataset.RasterCount, loc_dataset.RasterYSize, loc_dataset.RasterXSize))
    for line in range(loc_dataset.RasterYSize):
        # Read line in
        loc_data[:,line:line+1,:] = loc_dataset.ReadAsArray(0, line, loc_dataset.RasterXSize, 1)

    valid = np.logical_not(np.any(loc_data == nodata_value,axis=0))
    if trim_lines != 0:
        valid[:trim_lines, :] = False
        valid[-trim_lines:, :] = False

    # Grab zensor position and orientation information
    mean_latitude = find_angular_seeds(loc_data[1,valid].flatten(),1)
    mean_longitude = find_angular_seeds(-1 * loc_data[0,valid].flatten(),1)

    mean_elevation_km = np.mean(loc_data[2,valid]) / 1000.0

    # make elevation grid
    if lut_params.num_elev_lut_elements == 1:
        elevation_lut_grid = None
    else:
        min_elev = np.min(loc_data[2, valid])/1000.
        max_elev = np.max(loc_data[2, valid])/1000.
        elevation_lut_grid = np.linspace(max(min_elev, EPS),
                                         max_elev,
                                         lut_params.num_elev_lut_elements)

    return mean_latitude, mean_longitude, mean_elevation_km, elevation_lut_grid