I've been trying to find the best way to calculate the rank of a value from a 3d numpy array. This array is created from 35 years worth of rainfall data rasters. I'm treating the last raster in the stack as the "base" raster for comparison of the ranking in this case.

The numpy array's shape would be something like (36, 500, 500). For each column I need to get all the values from the stack and then rank where the most recent (last) year falls.

So far I have:

for row in range(rows):
    for col in range(cols):

        # get an array of all the values of the raster 
        # stack at current row and col iteration
        pixel_stack = stor_array[:, row, col]

        # Get the rank index of the last value in the pixel_stack
        rank_index = int(stats.rankdata(pixel_stack, 'min')[-1])

        # add the rank index at current 
        # row and col iteration to and output array
        out_array[row, col] = rank_index

This method works but seems slow. Is there a better (quicker/cleaner) way to iterate over the columns of the array?

For instance, is there anything that would allow me to use scipy function in a similar way to this numpy one?

out_arr = np.mean(stor_arr, axis=0)
  • Although numpy and scipy are used with Python alongside spatial libraries, this question seems to be using them in a way that has no spatial component so I think it would be better researched/asked at Stack Overflow. – PolyGeo Apr 29 '17 at 12:13

This function will do the ranking across an axis: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.mstats.rankdata.html

I'm not sure if your axis is 1 or 0, but give that a go.

(I don't know if there's one that will do that without the masking, but you can ignore that.)

| improve this answer | |
  • Thanks. I'll check that out at work tomorrow. Any idea if I could apply my own function in a similar manner if I wanted to calculate something else? – spalka Apr 28 '17 at 0:58
  • Sorry, I'm not sure how these functions work on the inside! I just use them =) – Alex Leith Apr 28 '17 at 1:21

After some more research I think I found the best solution for my needs. Numpy has an apply_along_axis function that let's you apply a function to a 1d slice of an input array. I was able to use this in conjunction with a lambda function and Scipy's rankdata function to achieve my desired result.

My resulting code is:

out_arr = np.apply_along_axis(lambda a: int(stats.rankdata(a, 'min')[-1]), 0, arr)

Hopefully someone finds this useful.

| improve this answer | |

Not the answer you're looking for? Browse other questions tagged or ask your own question.