如何在Python中找到二维数组的临界点？

Question

我有一个大小为(960,960)的数组，我想找出关键点，这样就能找到局部极值。

我尝试过使用np.diff和np.gradient这两个函数，但遇到了一些问题，不太确定该用哪个函数。

np.diff可以计算二阶差分，但gradient函数就没有这个选项。

我该怎么做才能找到关键点呢？

我试过

diff = np.diff(storm, n=2)                    

dxx = diff[0]                                                                                                                                  
dyy = diff[1]                                                                                                                                  

derivative = dyy/dxx 

在这里我遇到问题，因为dxx中的一些值等于零。

然后还有一个选项是

gradient = np.gradient(storm)
g2 = np.gradient(gradient)

但这个能帮我找到我想要的结果吗？

Answer 1

临界点就是一个函数的变化率（在多维情况下是梯度）等于0的地方。因此，你需要检查一下函数在x和y方向上的变化情况。numpy的diff函数在这种情况下非常好用。

如果在x和y方向上，相邻两个点的变化差值接近于0，那么你就可以说这个点是一个临界点。也就是说，当变化的方向发生变化（从负变正，或者从正变负）时，假设你的函数是平滑的。

# get difference in x- and y- direction
sec_grad_x = np.diff(storm,n=1,axis=0)
sec_grad_y = np.diff(storm,n=1,axis=1)

cp = []
# starts from 1 because diff function gives a forward difference
for i in range(1,n-1):
    for j in range(1,n-1):
        # check when the difference changes its sign
        if ((sec_grad_x[i-1,j]<0) != (sec_grad_x[i-1+1,j]<0)) and \
           ((sec_grad_y[i,j-1]<0) != (sec_grad_y[i,j-1+1]<0)):
            cp.append([i,j,  storm[i,j]])

cp = np.array(cp)