CF84 - 2 months ago 17x
Python Question

# Python: correct location of peaks in a 2D numpy array?

I know this question falls within the category of peak detection and there are answers available, but I think my problem is pretty simplistic and refers to a proof of principle.

Say I generate multiple,

`nxn`
2D numpy arrays of float values like this, which refer to a regular grid of
`nxn`
points (discrete domain):

``````myarray=
array([[ 0.82760819,  0.82113999,  0.811576  ,  0.80308705,  0.81231903,
0.82296263,  0.78448964,  0.79308308,  0.82160627,  0.83475755,
0.8580934 ,  0.8857617 ,  0.89901092,  0.92479025,  0.91840606,
0.91029942,  0.88523943,  0.85798491,  0.84190422,  0.83350783,
0.83520675],
[ 0.84971526,  0.84759644,  0.81429419,  0.79936736,  0.81750327,
0.81874686,  0.80666801,  0.82297348,  0.84980788,  0.85698662,
0.87819988,  0.89572185,  0.89009723,  0.90347858,  0.89703473,
0.90092666,  0.88362073,  0.86711197,  0.84791422,  0.83632138,
0.83685225], ...] #you can generate any nxn array of random values
``````

Now let's normalize them:

``````peak_value=myarray.max()
norm_array=myarray/peak_value
``````

And proceed to locate the
`(x,y)`
coordinates of the peak:

``````from collections import Counter
longit=[] #x values of the peak
latit=[] #y values of the peak
for x in range(myarray.shape[0]):
for y in range(myarray.shape[1]):
if norm_array[x][y]==1:
longit.append(x)
latit.append(y)

x=numpy.array(longit)
y=numpy.array(latit)

c=zip(x,y)
temp=Counter(elem for elem in c) #Counts the number of peaks in each (x,y) point in the 11x11 grid
d=dict(Counter(temp)) #Has the shape d={(x,y): number of peaks}
``````

Now this is just a single realization of the 2D numpy array. Given multiple arrays, the question is:

Is this the correct way to find the (x,y) of the peaks? Is there a more efficient way? Consider that there might be multiple peaks.

In C/C++ it is considered dangerous to do `==` with floating point numbers.

If there are no multiple exactly same pikes, you can use `numpy.argmax`:

``````a = random.rand(13, 11)
idx_unrolled = argmax(a)
the_peak = a[idx_unrolled/11, idx_unrolled%11]
``````

If you need all pikes, you can get list of `i`, `j` indices using `numpy.where`:

``````i, j = where(a > 0.99*the_peak)
``````

Use required number of `9` to tune the margin. For single precision floating points it is not than close to `1`.

The best way may be something like [http://stackoverflow.com/a/19141711/774971 ]:

``````i, j = where(a > (1.0 - 5*np.finfo(a.dtype).eps)*the_peak)
``````