1

I have two np array,'array 1' and 'array2' , and I have one image with 3 bands.

I want to stack array 1 and array two to be one image with two bands with the same shape as my image and then to clip it with another shapefile that I have. Array 1 and Array 2 are coming originally from the image so there are enough pixels to be the same shape.

This is how I have tried to do that:

#Load original image and check its' shape:
img=rasterio.open("original.tif")
array=img.read()
array.shape

>>>(3,2199,4041)

Here I did many calculations that in the end created one pandas df when each pixel from Original image is a row withour the shape of original image. I have created two numpy arrays from this final table which I want to construct into one image with to bands in the end:

Array_one=df_matrix['one'].values.reshape(2199, 4041)
Array_two=df_matrix['two'].values.reshape(2199, 4041)

Then I stack the two arrays together:

stack=np.stack((Array_one,Array_two))

stack
>>>array([[[0.6, 0.7, 0.7, ..., 0.2, 0.2, 0.2],
        [0.6, 0.7, 0.7, ..., 0.2, 0.2, 0.2],
        [0.6, 0.7, 0.7, ..., 0.2, 0.2, 0.2],
        ...,
        [0.5, 0.5, 0.5, ..., 0.2, 0.2, 0.3],
        [0.5, 0.5, 0.5, ..., 0.2, 0.2, 0.3],
        [0.4, 0.4, 0.5, ..., 0.3, 0.2, 0.3]],

       [[0.3, 0.3, 0.3, ..., 0.3, 0.4, 0.3],
        [0.3, 0.2, 0.2, ..., 0.4, 0.4, 0.3],
        [0.2, 0.2, 0.2, ..., 0.3, 0.3, 0.2],
        ...,
        [0.3, 0.5, 0.5, ..., 0.2, 0.2, 0.2],
        [0.4, 0.4, 0.5, ..., 0.2, 0.2, 0.2],
        [0.3, 0.4, 0.4, ..., 0.2, 0.2, 0.2]]])

Now things start to get complicated. I wanted my numpy to be saved as raster so i'll be able to clip it with rasterio or gdal but I couldn't find a way to do it. also to visualize it failed all the time:


# Use bilinear interpolation (or it is displayed as bicubic by default).
imshow(stack, interpolation="nearest")  
show()

TypeError: Invalid shape (2, 2199, 4041) for image data

that also happens if I use show(1). so then I tried to save it as raster with two band so I can load it again and then clip it , kike this:

 with rasterio.open('test', 
                    'w',
                    driver='GTiff',
                    height=stack.shape[0],
                    width=stack.shape[1],
                    count=2,
                    dtype=stack.dtype,
                    crs=img.crs,
                    nodata=None, # change if data has nodata value
                    transform=img.transform) as dst:
        dst.write(stack, 1)

but it fails all the time with this error:

ValueError: Source shape (1, 2, 2199, 4041) is inconsistent with given indexes 1

I pretty much sure that there is something i'm missing regard the number of bands.

My endgoal is to create raster which has two bands- array1 and array 2, so I can clip it with my shapefile.

Edit: I have tried the solution suggested by urban87 , I could display the numpy array : enter image description here

but when I saved it as raster and load it again like this:

 with rasterio.open('saveThis.tif', 
                    'w',
                    driver='GTiff',
                    height=stack.shape[0],
                    width=stack.shape[1],
                    count=2,
                    dtype=stack.dtype,
                    crs=img.crs,
                    nodata=None, # change if data has nodata value
                    transform=img.transform) as dst:
        dst.write(stack[0], 1)
        dst.write(stack[1], 2)


img=rasterio.open("saveThis.tif")
show(img,1)

I got this raster:

enter image description here

2
show(stack[0])
show(stack[1])

will work fine.

If you choose tuple, you have to choose the channel to output. Multi-band output is possible only when the raster dataset is loaded.

The same is true for the creation of raster files.

dst.write(stack, 1) 

will change

dst.write(stack[0], 1)
dst.write(stack[1], 2)

UPDATE
The cause of the error is that the width and height of the image are incorrectly specified.

stack.shape

return (band, height, width).

so, it will change

height=stack.shape[1],
width=stack.shape[2]

You can also use this.

height=img.height,
width=img.width
| improve this answer | |

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.