Summing all rasters in folder using Python

Question

I have a folder containing many raster files (.tif). I want to sum all of the rasters in this folder into a single output raster using Python in Jupyter Notebook. I have been able to use GDAL command line commands in OSGeo4W Shell, but I am not sure how to format this command into Python code that I can use in Jupyter Notebook.

To begin with, here is the code I am using to produce my rasters:

polygons = gpd.read_file('Boroughs_Test/Boroughs.shp')
polygon_IDs = polygons['ID'].tolist()

for i in polygon_IDs:
    x = polygons.loc[polygons['ID'] == i]
    vector_fn = x
    out_grid = make_geocube(
        vector_data=vector_fn,
        measurements=["test_value"],
        resolution=(-25, 25),
        fill=-9999,
    )
    out_grid["test_value"].rio.to_raster(str(i) + "_Output_Raster.tif")

This above code is taking my Boroughs.shp polygons, each polygon representing each of the 5 NYC boroughs, and rasterizing each polygon and sending it to its own unique raster. This means that for one raster file you would just see a rasterized polygon of Manhattan, and for the next raster file you would just see a rasterized polygon for Brooklyn, etc. I assigned a made-up value to each Borough called "test_value", which is the value I want to sum across the rasters.

So far what I have for a GDAL command to sum the rasters is:

for %f in (*.tif) do gdal_calc -A %f -B Result_Raster.tif --outfile=Result_Raster.tif --calc=A+B

where "Result_Raster.tif" is just a blank raster file with all 0 values.

How could I change the syntax of my command so that it can work with the gdal_calc.py function, but usable in Jupyter Notebook?

I am thinking this might involve a Python API, but I am not sure what the syntax would look like here. I am open to using rasterio as well.

Answer 1

This should calculate the sum of a lot of rasters as long as they are all the same resolution, extent, CRS, etc. I put in an assert statement to double check that.

import glob
import numpy as np
import rasterio


all_files = glob.glob('./raster_folder/*.tif')

# Create an initial array
with rasterio.open(all_files[0]) as src:
    result_array = src.read()
    result_profile = src.profile 

# Add on the rest one at a time
for f in all_files[1:]:
    with rasterio.open(f) as src:
        # Only sum the arrays if the profiles match. 
        assert result_profile == src.profile, 'stopping, file {} and  {} do not have matching profiles'.format(all_files[0], f)
        result_array = result_array + src.read()
        
with rasterio.open('Result_raster.tif', 'w', **result_profile) as dst:
        dst.write(result_array, indexes=[1])