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.