1

I have 2 multiband rasters and I want to mask some areas of one of them with values of the other one, so I write this code (I can't figure out why is not working):

import os, rasterio
import numpy as np

feb = r'path\to\multiband\stack_1.img'
may = r'path\to\multiband\stack_2.img'
mask = r'path\to\mask_raster' #this one is just one band
outfile = r'path\to\outfile' #multiband

with rasterio.open(mask) as mascara:

    MASK = mascara.read() 

with rasterio.open(feb) as febrero:

    FEB = febrero.read([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) 


with rasterio.open(may) as mayo:

    MAY = mayo.read([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) 

for i in range(1,11):

    MAY = np.where(MASK == 1, FEB[i], MAY[i])

profile = mayo.meta
profile.update(dtype=rasterio.int16)

with rasterio.open(outfile, 'w', **profile) as dst:
    dst.write(MAY.astype(rasterio.int16))

I get this error when running the code:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-21-054d59b9e8d1> in <module>()
     22 for i in range(1,11):
     23 
---> 24     MAY = np.where(MASK == 1, FEB[i], MAY[i])
     25 
     26 profile = mayo.meta

IndexError: index 2 is out of bounds for axis 0 with size 1

I am able to do it one by one band, but then I have to do the stack again and overall I would like to know why the code doesn't works.

3

It looks to me that you need MAY[i] = np.where(MASK == 1, FEB[i], MAY[i]) in your code. Note the indexing on the left hand side. Without this, MAY[i] will have the wrong shape on your second (i=2) loop iteration.

Rasterio tip: get the base profile for writing your output file from profile = mayo.profile instead of meta, it's more future-proof.

  • Yes, it was the MAY[i] in the left side. 1000 thanks, I hope I am gonna save a lot of time now! Only one thing about the tip, when I change mayo.meta with mayo.profile it returns an error: "ValueError: can't read closed raster file", so I did it in the old fashion way! – Digd Jul 7 '16 at 18:49

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.