Creating a raster based on values of other rasters

Question

I have three single band rasters, raster_a, raster_b and raster_c with pixel values 1 and 0.

I am trying to burn the values of raster_a to raster_b and create a "new_raster_b" where if: raster_a == 1, raster_b becomes 1 and if raster_a ==0 raster_b becomes 0. Special condition: if raster_c==0 raster_b remains 0 even if raster_a ==1.

I am reading all the 3 rasters as gdal datasets and converting them to array. How can I apply the above mentioned conditions to create an array for new_raster_b using these arrays?

raster_a_ds = gdal.Open(raster_a_path) 
raster_b_ds = gdal.Open(raster_b_path) 
raster_c_ds = gdal.Open(raster_c_path) 

print('Reading all rasters as arrays')
dcl_array = raster_a_ds.ReadAsArray()
tcm_array = raster_b_ds.ReadAsArray()
tcd_array = raster_c_ds.ReadAsArray()

Answer 1

Assuming your three rasters have the same dimensions, you can use numpy's boolean indexing to accomplish this. First, you need to create three masks, each one corresponding to one of your conditions:

con1 = (dcl_array == 1)  # raster a is 1
con2 = (dcl_array == 0)  # raster a is 0
con3 = (tcd_array == 0)  # raster c is 0

Then, you just have to index the original array using these masks and assign the values you mentioned:

tcm_array[con1] = 1  # replace with 1 where a is 1
tcm_array[con2] = 0  # replace with 0 where a is 0
tcm_array[con3] = 0  # replace with 0 where c is 0

Note that your special condition is executed last in order to make sure those cells remain 0.

Finally, you just have to get the raster band and write the modified array to it.

band = raster_b_ds.GetRasterBand(1)
band.WriteArray(tcm_Array)
band.FlushCache()

del band, raster_a_ds, raster_b_ds, raster_c_ds

Note that you might need to open the raster in writing mode if you want to edit its content. Simply change your second line to:

raster_b_ds = gdal.Open(raster_b_path, 1)  # 1 is for opening in writing mode

As I don't have access to your data and did not test this, I strongly recommend you to try this code on a copy before overwriting your original raster.

---

Regarding your comment, I'd copy the original array and then compare it to the modified one so you can get a mask to select all the modified values. Then, using the same logic of above, I'd assing the new value to it.

import numpy as np
# ...
# all your code from above
tcm_array = raster_b_ds.ReadAsArray()
tcm_array_original = np.copy(tcm_array)

# the code I wrote above
con1 = (dcl_array == 1)  # raster a is 1
# ...
# ...
# del band, raster_a_ds, raster_b_ds, raster_c_ds

# new logic
diff = (tcm_array != tcm_array_original)
tcm_array_original[diff] = 102

# save to the new raster
ds = gdal.Open("path", 1)
band = ds.GetRasterBand(1)
band.WriteArray(tcm_array_original)

del ds, band