[GIS] Stack two numpy array into one image with two bands in order to clip it with anotherr shapefile

clipnumpypythonrasterrasterio

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

Best Answer

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
Related Question