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()

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)

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)

del ds, band
  • thanks for the solution it worked. I am stuck at the next step now in I have to burn the changed values of tcm_array to another raster_p which has values from 0-101. The changed values of tcm_array has to be burned as value of 102. Can you please help me with building the logic for that? – FJ_Abbasi Jan 19 at 18:48
  • Sure, I just edited my answer. Take a look and let me know if it works. – Marcelo Villa Jan 19 at 19:17
  • Thanks a lot. It works :) – FJ_Abbasi Jan 19 at 20:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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