1

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.

2
  • Are all your raster files the same resolution and CRS?
    – Shawn
    Dec 17, 2021 at 1:50
  • Yes, all within the same CRS since they are all rasterized polygons produced from the same shapefile. I am trying to figure out how to loop through the folder and take raster_1, and then add raster_2 to that, and then add raster_3 to that, etc. I am having trouble figuring out simple rasterio code that can accomplish this. Dec 17, 2021 at 1:53

1 Answer 1

2

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])
5
  • Thanks, this looks promising, though I tried it out and received this error: AssertionError: stopping, file ./raster_folder\Raster_1.tif and ./raster_folder\Raster_2.tif do not have matching profiles. I am hoping this doesn't mean they have a different CRS, I don't see why they would, but I don't know. The traceback points to this line as the issue: assert result_profile == src.profile, 'stopping, file {} and {} do not have matching profiles'.format(all_files[0], f) Dec 17, 2021 at 2:14
  • It means there is something different about them that makes a straightfoward addition not so easy. Can you do this for each to see the difference: print(rasterio.open('./raster_folder\Raster_1.tif').profile)
    – Shawn
    Dec 17, 2021 at 2:19
  • Ok, here are the print results. Raster_1.tif: {'driver': 'GTiff', 'dtype': 'float64', 'nodata': -9999.0, 'width': 467, 'height': 877, 'count': 1, 'crs': CRS.from_epsg(32618), 'transform': Affine(25.0, 0.0, 580450.0, 0.0, -25.0, 4525925.0), 'tiled': False, 'interleave': 'band'}, Raster_2.tif: {'driver': 'GTiff', 'dtype': 'float64', 'nodata': -9999.0, 'width': 567, 'height': 578, 'count': 1, 'crs': CRS.from_epsg(32618), 'transform': Affine(25.0, 0.0, 589900.0, 0.0, -25.0, 4529975.0), 'tiled': False, 'interleave': 'band'} Dec 17, 2021 at 2:24
  • The crs is the same but the transform and width/height are not. So they are not exactly aligned. If you looked at them in ArcGIS or QGIS they would be slightly overlapping. Either here or in the rasterize from polygon step you need to reproject each to a common size and extent like in this answer stackoverflow.com/a/59109386
    – Shawn
    Dec 17, 2021 at 13:58
  • I had not considered this, but thanks for pointing me to this post. For context, each raster file represents a rasterized polygon of each of the 5 NYC boroughs. I updated my post with my initial rasterizing code to show how I generated these rasters I am trying to sum. I am just unsure why they would not be aligned since they were all produced from the same NYC shapefile. I looked at the rasters in QGIS, and I am not sure about overlapping, since they just fit together where they should be to make up NYC like puzzle pieces. I'm wondering if I will need to fix my initial rasterizing code. Dec 17, 2021 at 17:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.