QGIS Python Script – Using Raster Calculator in QGIS Python Script

bandpyqgisqgisrasterraster-calculator

I have a large set of radar .tif files with several bands, and I would like to have some Python scripts to perform some basic operations that I would usually perform using the Raster Calculator in QGIS 3.22.6, namely:

Create a new set of .tif files with only one of the bands of the input files, with the name equal to the original raster with a suffix (like here)
Create a new set of .tif files in which the output is a single band resulting from applying (a+b+c+d)/4 to the first 4 bands of the input files.
Create a new set of .tif files in which the output is a single band resulting from applying d/b to the 2nd and 4th input file bands.

I am currently trying to do this based on this similar question, but I don't know how to use other layers than the first one. My attempt did not produce the expected results:

from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator
from qgis.core import QgsProject

lddLrs = [tree_layer.layer() for tree_layer in QgsProject.instance().layerTreeRoot().findLayers()]

path = "path/to/output"

for lyr in lddLrs:
     entries = []
     ras = QgsRasterCalculatorEntry()
     ras.ref = 'ras@1'
     ras.raster = lyr
     ras.bandNumber = 1
     entries.append(ras)
     calc = QgsRasterCalculator('("infra-red@1"+"infra-red@2"+"infra-red@3"+"infra-red@4")/4', path + lyr.name() + "_suffix.tif", 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries)

calc.processCalculation()

This other question seemed to have a similar problem, but there doesn't seem to be an answer to it yet.

Best Answer

I think it's easier to use GDAL raster calculator. I execute the tool manually once in QGIS, press Ctrl+Alt+H and copy the command. Then write the loop.

To create the output names you can use layer.source() and os.path.join with os.path.basename.

For first band, and mean:

import os
output_folder = r'/home/bera/Desktop/GIStest/rasters/'
for layer in QgsProject.instance().mapLayers().values(): #For all layers in the map
    print(layer)

    #########Extract first band
    #Create the output name
    #layer.source() 
        #'/home/bera/Desktop/GIStest/4bandraster.tif'
    #os.path.basename(layer.source()) #
        #'4bandraster.tif'
    output1 = os.path.join(output_folder, os.path.basename(layer.source()).replace('.tif', '_firstbandonly.tif'))
        #'/home/bera/Desktop/GIStest/rasters/4bandraster_suffix.tif'
    processing.run("gdal:rastercalculator", {'INPUT_A':layer,'BAND_A':1,'INPUT_B':None,
        'BAND_B':None,'INPUT_C':None,'BAND_C':None,'INPUT_D':None,'BAND_D':None,
        'INPUT_E':None,'BAND_E':None,'INPUT_F':None,'BAND_F':None,'FORMULA':'A',
        'NO_DATA':None,'RTYPE':5,'OPTIONS':'','EXTRA':'','OUTPUT':output1})
    
    ########Mean of bands 1,2,3,4:
    output2 = os.path.join(output_folder, os.path.basename(layer.source()).replace('.tif', '_mean1234.tif'))
    processing.run("gdal:rastercalculator", {'INPUT_A':layer,'BAND_A':1,'INPUT_B':layer,'BAND_B':2,'INPUT_C':layer,'BAND_C':3,'INPUT_D':layer,'BAND_D':4,'INPUT_E':None,'BAND_E':None,'INPUT_F':None,'BAND_F':None,'FORMULA':'(A+B+C+D)/4','NO_DATA':None,'RTYPE':5,'OPTIONS':'','EXTRA':'','OUTPUT':output2})
    
print("Done")

Related Solutions

GDAL – How to Use gdal_retile on RGBI (4 Band) GeoTIFFs While Preserving Band Information

I found an answer for exactly this problem in the [gdal-dev] Mailing list: To avoid the fourth band being interpreted as ALPHA channel, the corresponding option ALPHA=NO can be specified for the GeoTIFF driver. To preserve the RGB order, apparently there is the option PHOTOMETRIC=RGB - for more read the mailing list. In my case the ALPHA=NO option was enough.

When calling the script via command-line, use the -co "OPTION-NAME=VALUE" parameter (see doc) to specify options.

My full command was:

gdal_retile.py 

-co "ALPHA=NO"

-s_srs EPSG:31464 [georeference]
-ps 1024 1024 [pixel width/height)
-tileIndex "tileIndex.shp" 
-csv "tileGeoRef.csv" 
-targetDir "...\outputfolder" 
"...\inputfile.tif"

Python – Open Sentinel 2 Imagery from AWS Servers Using Rasterio Module

You are reading and thus downloading the entire image for both of those bands (B08.tif and B04.tif).

It is much quicker to only read the subset of those images that overlap with the polygon instead of clipping it out later on after downloading the lot. This is one of the primary use cases for the Cloud-Optimized Geotiff (COG).

The below example uses rasterio windowed reading and took 6.5 seconds total for all 3 polygons in my test dataset.

import rasterio
import geopandas as gpd
from rasterio.windows import transform
from rasterio.features import geometry_window, geometry_mask

# Sentinel-2 bands
nir_path = 's3://sentinel-cogs/sentinel-s2-l2a-cogs/36/N/YF/2023/6/S2B_36NYF_20230605_0_L2A/B08.tif'
red_path = 's3://sentinel-cogs/sentinel-s2-l2a-cogs/36/N/YF/2023/6/S2B_36NYF_20230605_0_L2A/B04.tif'

# Polygon path (the file contains 3 polygons with a unique ID field)
polygon_path = 'path_to_sample_polygon.geojson'
polygon_id_field = "ID"

# Load and transform the polygon
polygon_gdf = gpd.read_file(polygon_path)
polygon_gdf = polygon_gdf.to_crs("EPSG:32636")

# Open the AWS COG datasets
with rasterio.Env(AWS_NO_SIGN_REQUEST="YES"):
    with rasterio.open(nir_path) as nir_src, rasterio.open(red_path) as red_src:
        
        # loop through each feature in the GDF
        for feature in polygon_gdf.iterfeatures(show_bbox=True):

            # Get a window
            window = geometry_window(nir_src, [feature["geometry"]])
            window_transform = transform(window, nir_src.transform)
            window_shape = (window.height, window.width)

            # Read all the data in the window, masking out any NoData
            nir = nir_src.read(window=window, masked=True).astype('float32')
            red = red_src.read(window=window, masked=True).astype('float32')

            # Update the NoData mask to exclude anything outside the polygon
            mask = geometry_mask([feature["geometry"]], window_shape, window_transform)
            nir.mask += mask
            red.mask += mask

            # Calculate NDVI
            ndvi = (nir - red) / (nir + red)

            # Save the masked NDVI to one tif per polygon
            ndvi_path = f'ndvi_{feature["properties"][polygon_id_field]}.tif'

            meta = nir_src.meta
            # Don't forget to update the dtype! Your original code didn't, so the Int16 output truncated NDVI values to all 0's
            meta.update({"driver": "GTiff", "dtype": ndvi.dtype, "height": window.height, "width":window.width, "transform": window_transform})
            with rasterio.open(ndvi_path, 'w', **meta) as dst:
                dst.write(ndvi)

Best Answer

Related Solutions

GDAL – How to Use gdal_retile on RGBI (4 Band) GeoTIFFs While Preserving Band Information

Python – Open Sentinel 2 Imagery from AWS Servers Using Rasterio Module

Related Question