How to use the gdal.GA_Update function in GDAL

import os.path
try:
    from osgeo import gdal
except ImportError:
    import gdal
#@-<>
#@+others
#@+node:maphew.20100601093031.2401: ** usage
if len(sys.argv) < 2:
    print "Usage: gdalsetnull.py raster_file null_value {null band2} {null band3} ..."
    sys.exit(1)
#@+node:maphew.20100601093031.2400: ** Main
input = sys.argv[1]
null_value = sys.argv[2]

dataset = gdal.Open( input, gdal.GA_Update )
if dataset is None:
    print 'Unable to open', input, 'for writing'
    sys.exit(1)

for band_num in range(1, dataset.RasterCount+1):
    band = dataset.GetRasterBand(band_num)
    print 'Initial nodata for band %s \t %s'  % (band_num,band.GetNoDataValue() )

    # the 2nd commandline argument is for band#1, arg3 for band2, etc.
    arg_num = band_num + 1 
    # optionally set different nodata values for each bands
    if arg_num > 2:
        if sys.argv[arg_num]:
            null_value = sys.argv[band_num]
            
    # FIXME: handle case where we want to remove nodata altogether (change to 'None')

def transfer_metadata(aligned_image_path, original_image_path):
    aligned_dataset = gdal.Open(aligned_image_path, gdal.GA_Update)
    original_dataset = gdal.Open(original_image_path, gdal.GA_ReadOnly)

    aligned_dataset.SetMetadata(original_dataset.GetMetadata())
    rpcs = original_dataset.GetMetadata('RPC')
    aligned_dataset.SetMetadata(rpcs, 'RPC')
    aligned_dataset.SetGeoTransform(original_dataset.GetGeoTransform())
    aligned_dataset.SetProjection(original_dataset.GetProjection())

    aligned_dataset = None
    original_dataset = None

considered valid.

    Args:
        infile: string or os.path object for a
                gdal compatible intensity file of Sentinel-1

    Notes:
        The file will be manipulated inplace, meaning,
        no new outfile is created.
    '''

    # print(' INFO: Removing the GRD Border Noise.')
    currtime = time.time()

    # read raster file and get number of columns adn rows
    raster = gdal.Open(infile, gdal.GA_Update)
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    # create 3000xrows array for the left part of the image
    array_left = np.array(raster.GetRasterBand(1).ReadAsArray(0,
                                                              0, 3000, rows))

    for x in range(3000):
        # condition if more than 50 pixels within the line have values
        # less than 500, delete the line
        # if np.sum(np.where((array_left[:,x] < 200)
        # & (array_left[:,x] > 0) , 1, 0)) <= 50:
        if np.mean(array_left[:, x]) <= 100:
            array_left[:, x].fill(0)
        else:
            z = x + 150

raster_class_name = f"{raster_name}_class.tif"
        raster_class_path = os.path.join(output_path, raster_class_name)
        # Create the class raster output
        class_ds = gdal.GetDriverByName('GTiff').Create(raster_class_path, xsize, ysize, 1, gdal.GDT_Byte)
        if class_ds is None:
            raise IOError(f"Unable to create: [{raster_class_path}]")
        else:
            logger.debug(f"Creating: [{raster_class_path}]")
        class_ds.SetGeoTransform(affine)
        class_ds.SetProjection(georef)
        class_band = class_ds.GetRasterBand(1)
        class_band.SetNoDataValue(0)
        class_ds.FlushCache()   # save to disk
        class_ds = None     # noqa # need to close before open-update
        class_band = None   # noqa # also close band (remove ptr)
        class_ds = gdal.Open(raster_class_path, gdal.GA_Update)
        # Create the probabilities raster output
        raster_prob_name = f"{raster_name}_probs.tif"
        raster_prob_path = os.path.join(output_path, raster_prob_name)
        probs_ds = gdal.GetDriverByName('GTiff').Create(raster_prob_path, xsize, ysize, class_count, gdal.GDT_Float32)
        if probs_ds is None:
            raise IOError(f"Unable to create: [{raster_prob_path}]")
        else:
            logger.debug(f"Creating: [{raster_prob_path}]")
        probs_ds.SetGeoTransform(affine)
        probs_ds.SetProjection(georef)
        probs_ds.FlushCache()   # save to disk
        probs_ds = None  # noqa # need to close before open-update
        probs_ds = gdal.Open(raster_prob_path, gdal.GA_Update)
        return class_ds, probs_ds

def buffer_mask_in_place(mask_path, buffer_size):
    """Expands a mask in-place, overwriting the previous mask"""
    log = logging.getLogger(__name__)
    log.info("Buffering {} with buffer size {}".format(mask_path, buffer_size))
    mask = gdal.Open(mask_path, gdal.GA_Update)
    mask_array = mask.GetVirtualMemArray(eAccess=gdal.GA_Update)
    cache = morph.binary_erosion(mask_array.squeeze(), selem=morph.disk(buffer_size))
    np.copyto(mask_array, cache)
    mask_array = None
    mask = None

def MeanMap():
    ''' 
    Create Arcus Map
    
    Treshhold the percentage per day map to 30%
    
    '''
    infile =  'C:\\Users\\max\\Documents\\Icecharts\\Arcus\\icechart_processed.tif'
    outfile = 'C:\\Users\\max\\Documents\\Icecharts\\Arcus\\arcus_map.tif'
    outshape = 'C:\\Users\\max\\Documents\\Icecharts\\Arcus\\arcus_map.shp'
    
    #Open Rasterfile and Mask
    driver = gdal.GetDriverByName('GTiff')
    driver.Register()
    ds = gdal.Open(infile, gdal.GA_Update)

    
    #Read input raster into array
    iceraster = ds.ReadAsArray()
    
    #get image max and min and calculate new range
    rows = ds.RasterYSize
    cols = ds.RasterXSize
    
    #create output images
    driver = ds.GetDriver()
    outraster = driver.Create(outfile, cols, rows, 1, gdal.GDT_Float64 )
    if outraster is None: 
        print 'Could not create ', outfile
        return

The routine checks the outer 3000 columns for its mean value.
    If the mean value is below 100, all values will be set to 0,
    otherwise the routine will stop, assuming all columns towards
    the inner image are valid.


    :param inFile: gdal compatible intensity file of Sentinel-1
    :return:
    '''

    print(' INFO: Removing the GRD Border Noise.')
    currtime = time.time()

    # read raster file and get number of columns adn rows
    raster = gdal.Open(inFile, gdal.GA_Update)
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    # create 3000xrows array for the left part of the image
    array_left = np.array(raster.GetRasterBand(1).ReadAsArray(0,
                                                              0, 3000, rows))

    for x in range(3000):
        # condition if more than 50 pixels within the line have values
        # less than 500, delete the line
        # if np.sum(np.where((array_left[:,x] < 200)
        # & (array_left[:,x] > 0) , 1, 0)) <= 50:
        if np.mean(array_left[:, x]) <= 100:
            array_left[:, x].fill(0)
        else:
            z = x + 150

def copy_tif_info(in_fpath, out_fpath):
    """
    Copy the TIF metadata from in_fpath and save it to the image at out_fpath.
    :param in_fpath: Filepath for input image which has metadata to copy.
    :type in_fpath: str

    :param out_fpath: Filepath for image to which we want to save the metadata.
    :type out_fpath: str
    """
    try:
        logging.info("---- Copying TIF information ----")
        gdal_in = gdal_utils.gdal_open(in_fpath)
        gdal_out = gdal_utils.gdal_open(out_fpath, gdal.GA_Update)
        metadata_domains = gdal_in.GetMetadataDomainList()
        for domain in metadata_domains:
            dico = gdal_in.GetMetadata_Dict(domain)
            for key, val in dico.items():
                gdal_out.SetMetadataItem(key, val, domain)
        return True
    except Exception as e:
        logging.exception(e)
        return False

directory=output_dir)

		print("Standardizing band " + str(band) + ' ' + image_path + " => " + output_image_path)

		if not Path(output_image_path).is_file():
			dataType = gdal.GDT_Float32
			nbands = len(bands)
			
			if convert_int16:
				dataType = gdal.GDT_Int16

			output_ds = dl_utils.create_output_file(image_path, output_image_path, \
				nbands, dataType)

		else:
			output_ds = gdal.Open(output_image_path, gdal.GA_Update)
		
		input_ds = gdal.Open(image_path, gdal.GA_ReadOnly)
		
		x_size = input_ds.RasterXSize
		y_Size = input_ds.RasterYSize

		for xoff in range(0,x_size,chunk_x_size):
			if (xoff+chunk_x_size) > x_size:
				chunk_x_size = x_size - xoff

			output_band_ds = output_ds.GetRasterBand(band)

			intput_band_ds = input_ds.GetRasterBand(band)
			band_data = intput_band_ds.ReadAsArray(xoff, 0, chunk_x_size, y_Size);
			band_data = band_data.astype('Float32')

dem_path
        The path to the DEM
    out_raster_path
        The path to the output.
    raster_datetime
        A datetime.DateTime object **with timezone set**

    """
   
    with TemporaryDirectory() as td:

        in_raster = gdal.Open(raster_path)
        in_array = in_raster.GetVirtualMemArray()
        out_raster = ras.create_matching_dataset(in_raster, out_raster_path, bands=in_raster.RasterCount,
                                                 datatype=gdal.GDT_Float32)
        out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)

        print("Preprocessing DEM")
        # Right, so.
        # Something in there is going funny, and I'm around 80% sure that it's because the DEM isn't aligning
        # properly with the
        clipped_dem_path = p.join(td, "clipped_dem.tif")
        reproj_dem_path = p.join(td, "reproj_dem.tif")
        ras.reproject_image(dem_path, reproj_dem_path, in_raster.GetProjection(), do_post_resample=False)
        ras.resample_image_in_place(reproj_dem_path, in_raster.GetGeoTransform()[1])  # Assuming square pixels
        ras.align_image_in_place(reproj_dem_path, raster_path)
        ras.clip_raster_to_intersection(reproj_dem_path, raster_path, clipped_dem_path, is_landsat)

        ic_array, zenith_array, slope_array = calculate_illumination_condition_array(clipped_dem_path, raster_datetime)
        
        if is_landsat:
            in_array = in_array.T