Search code examples
pythonnumpyrasterrasterio

Loop through each pixel of a raster stack and return a time series of pixels using Python


I'm new to Python language and I'm trying to loop through a rasterstack and store pixels in a time-series manner.

For example, suppose I have three rasters of three dates, for 2020, 2021, 2022:

 A = array([[[0.2, 0.3, 0.4, 0.5,
              0.6, 0.7, 0.8, 0.9]],
            [[1.0, 1.1, 1.2, 1.3,
              1.4, 1.5, 1.6, 1.7]],
            [[1.8, 1.9, 2.0, 2.1,
              2.2, 2.3, 2.4, 2.5]]])

I would like to create a new array with arrays whose elements are displayed like:

  B = array([[0.2, 1.0, 1.8],
            [0.3, 1.1, 1.9],
            [0.4, 1.2, 2.0],
            ...
            [0.9, 1.7, 2.5]])

i.e., [0.2, 1.0, 1.8] is formed by the first element (0.2) which was the first element of first array of A, the second element (1.0) is the first element of second array of A, the third element (1.8) is the first element of third array of A. Then for the next array [0.3, 1.1, 1.9], the first element (0.3) is the second element of first array of A. The second element (1.1) is the second element of second array of A and so on.

Is there any easy way to do this without a lot of loops? To get some data:

    data = np.random.random((3, 4, 4)) 
    stack = np.dstack(data) #just to change to (4,4,3), number of images for last

Solution

  • You can transpose the original array with the code below, I resize after the transpose to get rid of the extra dimension but depending on your actual data you might not need/want it.

    arr = np.array([
        [[0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]],
        [[1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7]],
        [[1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5]],
    ])
    tr = arr.transpose()
    #tr = array([
    #    [[0.2, 1. , 1.8]],
    #    [[0.3, 1.1, 1.9]],
    #    [[0.4, 1.2, 2. ]],
    #    [[0.5, 1.3, 2.1]],
    #    [[0.6, 1.4, 2.2]],
    #    [[0.7, 1.5, 2.3]],
    #    [[0.8, 1.6, 2.4]],
    #    [[0.9, 1.7, 2.5]],
    #])
    reshaped = tr.reshape((8,3))
    #reshaped = array([
    #    [0.2, 1. , 1.8],
    #    [0.3, 1.1, 1.9],
    #    [0.4, 1.2, 2. ],
    #    [0.5, 1.3, 2.1],
    #    [0.6, 1.4, 2.2],
    #    [0.7, 1.5, 2.3],
    #    [0.8, 1.6, 2.4],
    #    [0.9, 1.7, 2.5],
    #])