[GIS] Efficient way to multiply specific raster with specific variable python

gdalpythonraster

I have 4 separate raster bands of the same image and 4 different coefficients (one for each band) stored as variables. I would like to iterate over the raster folder, and multiply each band with its corresponded coefficient. I can recognize each band by its file name. Let's say I have in folder c:\bands these bands:

sn2_b1.tif
sn2_b2.tif
sn2_b3.tif
sn2_b3.tif

And let's say these are my coefficients:

coeff_b1 = 10
coeff_b2 = 20
Coeff_b3 = 30
Coeff_b4 = 40

My goal is iterate over the bands, multiply each one with its corresponded coefficient and save the result as a new raster. My solution so far is quite cumbersome:

import gdal
import os
src = r"c:\bands"
dst = r"c:\cal_bands"
for i in os.listdir(src):
    if i.endswith(".tif"):
        orig_tif = os.path.join(src,i)
        band_ds = gdal.Open(os.path.join(src,i))
        band_DN = band_ds.GetRasterBand(1).ReadAsArray()  
        geotransform = band_ds.GetGeoTransform()
        geoproj = band_ds.GetProjection()
        xsize = band_ds.RasterXSize
        ysize = band_ds.RasterYSize
        band_DN = band_DN.astype('f')
        format1 = "GTiff"
        driver = gdal.GetDriverByName(format1)
        dst_datatype = gdal.GDT_Float32
        if i.split("_")[1] == 'b1.tif':
            cal_band =  band_DN  * coeff_b1
        elif i.split("_")[1] == 'b2.tif':
            cal_band =  band_DN  * coeff_b2
        elif i.split("_")[1] == 'b3.tif':
            cal_band =  band_DN  * coeff_b3 
        elif i.split("_")[1] == 'b4.tif':
            cal_band =  band_DN  * coeff_b4
        new_name = "cal_" + i
        output_path = os.path.join(dst,new_name)
        dst_ds = driver.Create(output_path ,xsize,ysize,1,dst_datatype)
        dst_ds.GetRasterBand(1).WriteArray(cal_band)
        dst_ds.SetGeoTransform(geotransform)
        dst_ds.SetProjection(geoproj)
        dst_ds = None

This script works great but is there a way to make it less cumbersome? In this example I have only 4 bands, which is easy to write all the ifs, but I might have 20 bands or more.

Best Answer

It looks as though you are doing different operations inside a single loop, which you can break off into functions.

If I understand correctly, you:

Have a band-coefficient lookup table (LUT)
Unknown number of tifs in a source folder
Image band number is determined by the characters after the final underscore in the file name.

And you can match them together by extracting the band name from the string. The natural way to write LUTs in Python is with built-in dictionaries.

Code

import gdal
import os

src = r'c:\bands'
dst = r'c:\cal_bands'

COEFF_LUT = {'b1': 10, 'b2': 20, 'b3': 30, 'b4': 40}
FORMAT = 'GTiff'


def band_num(f):
    """Given raster filename with band number info, returns band number"""
    filename = f.split('.')[0]
    band = filename.split('_')[-1]
    return band


def multiply_raster_by_coeff(file, band):
    """Given filename and band number of image, multiplies it by saves a copy
    in specified destination folder"""
    coeff = COEFF_LUT[band]
    dst_datatype = gdal.GDT_Float32
    driver = gdal.GetDriverByName(FORMAT)

    # Read file, get metadata
    band_ds = gdal.Open(os.path.join(src, file))
    band_DN = band_ds.GetRasterBand(1).ReadAsArray()
    geotransform = band_ds.GetGeoTransform()
    geoproj = band_ds.GetProjection()
    xsize = band_ds.RasterXSize
    ysize = band_ds.RasterYSize

    band_DN = band_DN.astype('f')

    # Multiply by band coefficient
    cal_band = band_DN * coeff

    # Write to disk
    new_name = 'cal_' + file
    output_path = os.path.join(dst, new_name)
    dst_ds = driver.Create(output_path, xsize, ysize, 1, dst_datatype)
    dst_ds.GetRasterBand(1).WriteArray(cal_band)
    dst_ds.SetGeoTransform(geotransform)
    dst_ds.SetProjection(geoproj)


if __name__ == '__main__':
    files = [f for f in os.listdir(path) if f.endswith('.tif')]
    image_bands = {f: band_num(f) for f in files}
    for file, band in image_bands:
        multiply_raster_by_coeff(file, band)

Conclusion

You can now handle an arbitrary number of bands just by updating the COEFF_LUT dictionary.

You could further break the multiply_raster_by_coeff function. Of course, performance-wise the two scripts are equivalent, but future you (and collaborators) will appreciate this if it needs changes.

Related Solutions

[GIS] ny efficient way to create a raster with pixel value representing a unique id calculated based on pixel lat/lon

There are many ways to assign unique identifiers to raster cells. Some are listed below. Please note that although unique identifiers can be efficiently generated, any algorithm that requires unique identifiers for cells in large rasters (with many millions of more of cells) is itself going to be inefficient; in many cases it likely could be replaced by an algorithm not requiring every cell to have an identifier. Note, too, there is a natural limit to the size of such a grid, because even 32-bit grids can support at most 2^32-1 (about four billion) unique values. Methods to overcome that limit (without resorting to 64-bit grids) are given at the end.

Using a grid of column identifiers I (with values running from 0 up to one less than the number of columns, m) and a similar grid of row identifiers J (with values from 0 through n-1), compute
```
n * "I" + "J"
```
A convenient alternative is to round n up to the next power of ten equal to or exceeding n. The row and column coordinates can then be read directly off the point identifiers. For instance, with n = 12796 and m = 8092, use the formula 100000 * I + J. (Because the largest possible value of 809112795 will not overflow a signed 32-bit integer, this will work.) For instance, if a point has the identifier 49160008 it must have (row, column) coordinates (491, 60008) found by splitting the identifier into its last five digits and the remaining prefix thereof.

Column and row identifier grids can be computed (in one single calculation each) using the methods described at How do I calculate a weighted centroid with respect to a raster, e.g. population weighted centroid.
For smallish grids, simply use Combine("I", "J"). This operation will create internal unique pixel identifiers.
For very small grids (with fewer than 2^16 cells or so), generate a grid of uniformly distributed integers between 0 and the largest value possible (2^31-1 usually). Query the attribute table for counts larger than 1. If they exist, just generate a new random grid, repeating until all values are unique.
Don't even bother. If you ever need to identify a pixel, extract the values from the two grids I and J. Now there is no practical limit to the size of the grid--it need only have fewer than four billion rows and four billion columns.
The roles played by I and J can be replaced by grids of X and Y coordinates, respectively, but (a) the calculation in (1) is more delicate because it will be difficult to avoid potential collisions of IDs and (2) when the cellsize is very small (such as a few meters or less) single-precision floats usually won't distinguish all neighboring cells, so double precision will be necessary.
X and Y coordinate grids can themselves be replaced by other grids. Good ones are (signed) distances from appropriate curves: often two such distances suffice to uniquely identify every point close to the intersection of those curves. This provides a way to use non-Cartesian coordinate systems (such as polar or hyperbolic coordinates).

[GIS] Export raster layer in python to specific folder

Your code doesn't work because you're not doing anything in your for loop (except setting a file name).

Saving a raster is a bit tedious, afaik. Here is a sample code I gathered from varied sources; I didn't test it, so it might need some adjustments. It relies on the QgsRasterPipe and QgsRasterFileWriter classes.

layers = iface.layerTreeView().selectedLayers()
for layer in layers:
    file_name = '/full/raster/path/' + layer.name() + '.tif'
    file_writer = QgsRasterFileWriter(file_name)
    pipe = QgsRasterPipe()
    provider = layer.dataProvider()

    if not pipe.set(provider.clone()):
        print "Cannot set pipe provider"
        continue

    file_writer.writeRaster(
        pipe,
        provider.xSize(),
        provider.ySize(),
        provider.extent(),
        provider.crs())

Note: I'd love to know if there is a simpler way to save rasters!

EDIT: added dummy path in file_name (full path is required).