Convert 3D numpy to array to 4D array without changing

Question

I have a raster image with 10 bands and read it as array.

file = "test.tif"
ds = gdal.Open(file)
arr = ds.ReadAsArray()

The resulting shape of the 3D array looks like this: (n_bands, y_pixels, x_pixels)

However, the software I want to use requires a 4D array as input: (n_images, n_pixels_y, n_pixels_x, n_bands)

Is there a way to read the raster as array with the specified properties of a 4D array or to convert the 3D array to a 4D one.

I tried to use np.reshape, but it changes the location of pixels.

array([[[3344, 3344, 3344, ..., 8001, 8001, 8001],
    [3344, 3344, 3344, ..., 8001, 8001, 8001],
    [3344, 3344, 3344, ..., 8001, 8001, 8001],
    ...,
    [2359, 2359, 2359, ..., 7106, 7106, 7106],
    [2359, 2359, 2359, ..., 7106, 7106, 7106],
    [2359, 2359, 2359, ..., 7106, 7106, 7106]],
   ...,
    [[3173, 3173, 3431, ..., 5658, 5463, 5463],
    [3173, 3173, 3431, ..., 5658, 5463, 5463],
    [3393, 3393, 3487, ..., 5767, 5536, 5536],
    ...,
    [1751, 1751, 1722, ..., 2753, 2534, 2534],
    [1395, 1395, 1415, ..., 2672, 2521, 2521],
    [1395, 1395, 1415, ..., 2672, 2521, 2521]]], dtype=uint16)

arrn=arr.reshape(1,y_pixels,x_pixels,10)

array([[[[3344, 3344, 3344, ..., 2122, 2122, 2122],
     [2122, 2122, 1378, ..., 1378, 1420, 1420],
     [1420, 1420, 1420, ..., 1435, 1435, 1435],
     ...,
     [8753, 8753, 8753, ..., 8086, 8086, 8086],
     [8086, 8086, 6949, ..., 6949, 7091, 7091],
     [7091, 7091, 7091, ..., 7633, 7633, 7633]],
    ...,
     [[1944, 1944, 1885, ..., 1846, 1795, 1795],
     [1645, 1645, 1366, ..., 1605, 1706, 1706],
     [1723, 1723, 1854, ..., 2182, 2270, 2270],
     ...,
     [3057, 3057, 3059, ..., 3150, 3195, 3195],
     [3249, 3249, 3180, ..., 3178, 3165, 3165],
     [3145, 3145, 3056, ..., 2672, 2521, 2521]]]], dtype=uint16)

Answer 1

You want to move the n_bands axis to the end, and add a dimension in front. .reshape can't know that you want to do that and will just re-interpret the data in the new shape. But you can manually separate it into two steps to keep the correct order of your pixels:

arr  # shape (n_bands, y_pixels, x_pixels)
swapped = np.moveaxis(arr, 0, 2)  # shape (y_pixels, x_pixels, n_bands)
arr4d = np.expand_dims(swapped, 0)  # shape (1, y_pixels, x_pixels, n_bands)