Divide and write raster in GDAL using Python?

Question

I have a folder of raster .tif. files. I want to iterate through all the rasters in this folder and divide them by a specific raster. I then want to take these outputs and write them to new rasters.

Here is what I am trying:

d = ['Blue', 'Red', 'Green']
test = 'test.tif'
folder = Path('Path/to/My/Folder')

for f in folder.glob("*.tif"):
    if any(color in f.name for color in d):
        new_raster = f/test

What I am doing here is this: In my folder, I have many raster .tif files. However, I only actually want to use them in my loop code if their file names contain "Blue", "Red" or "Green". So I would want to use the files Blue_Raster.tif, Red_Raster.tif, and Green_Raster.tif in my loop, but exclude a Purple_Raster.tif, etc. So then FOR EACH of these rasters, I want to divide them by my selected raster test.tif. This produces FOR EACH of my rasters (just for Blue, Red, and Green) a new_raster, which I want to write to an output raster, titled either Blue_Output_raster.tif, Red_Output_raster.tif, or Green_Output_raster.tif. However, I am not sure how to actually write these rasters. Would this just involve a simple use of the writeraster() function in gdal? I am not sure if there are actually many more steps involved or if I even used the raster division here correctly. I am open to using rasterio or geocube as well if that is simpler.

Answer 1

This should work with what you require. Read about rasterio write here

d = ['Blue', 'Red', 'Green']
test = 'test.tif'
folder = Path('Path/to/My/Folder')

with rasterio.open("test.tif") as src:
  
  for f in folder.glob("*.tif"):
    if any(color in f.name for color in d):
      
      with rasterio.open(f) as med:    
        new_raster = med.read(1)/src.read(1)

      # write raster
      profile = src.profile
      profile.update(
          dtype=new_raster.dtype,
          count=1,
          compress='lzw')
      filename = "_output_".join(os.path.splitext(f)[0].split("_"))
      with rasterio.open(f'{filename}.tif', 'w', **profile) as dst:
          dst.write(new_raster, 1)

Answer 2

rioxarray can simplify that:

https://corteva.github.io/rioxarray/stable/examples/reproject_match.html

https://www.earthdatascience.org/courses/use-data-open-source-python/intro-raster-data-python/raster-data-processing/subtract-rasters-in-python/

https://carpentries-incubator.github.io/geospatial-python/07-raster-calculations/index.html

Answer 3

# Look at this simple code:

import os
import gdal


path = r'inputpath_directory'


# for Loop
 
bands = []

for i in os.listdir(path):
    bands.append(i)
    
print(bands) # a list of bands


# To save a raster using gdal:


gdal_array.SaveArray(bands, 
                     r'outputpath_name', 
                     "GTiff", 
                      prototype=ds) # band prototype