如何遍历numpy数组并选择相邻单元格

Question

我正在把美国地质调查局的高程栅格数据转换成一个Numpy数组，然后想随机选择数组中的一个位置。从这个位置开始，我想创建一个方法，找出周围八个单元格，看看这些单元格的高程是否在随机选择的单元格的一米范围内。

这就有点复杂了……如果一个邻居在一米之内，那么就会对它调用同样的方法，继续这个过程，直到没有单元格的高程在一米之内，或者选中的单元格数量达到了设定的限制。

如果这还不太清楚，希望下面这个二维数组的例子能让你更明白。加粗/斜体的单元格（35）是随机选择的，方法被应用到它上面（选择了它的八个邻居），然后这个方法又被应用到所有邻居上，直到没有更多的单元格可以选择（所有加粗的数字都是被选中的）。

33 33 33 37 38 37 43 40

33 33 33 38 38 38 44 40

36 36 36 36 38 39 44 41

35 36 35 35 34 30 40 41

36 36 35 35 34 30 30 41

38 38 35 35 34 30 30 41

我对Java还算不错，知道怎么写一个方法来实现这个目的，不过GIS主要是基于Python的。我正在学习Python，并且写了一些代码，但在把Python适应到GIS脚本接口时遇到了很大的问题。

谢谢你的帮助！

问题继续...

谢谢你的回答，Bas Swinckels。我试着把你的代码融入到我目前写的代码中，结果陷入了无限循环。下面是我写的代码。我需要克服两个主要步骤才能让它工作。这是我从栅格生成的数组的一个例子（-3.40e+38是无数据值）。

>>> 
[[ -3.40282306e+38  -3.40282306e+38  -3.40282306e+38 ...,  -3.40282306e+38
   -3.40282306e+38  -3.40282306e+38]
 [ -3.40282306e+38  -3.40282306e+38  -3.40282306e+38 ...,  -3.40282306e+38
   -3.40282306e+38  -3.40282306e+38]
 [ -3.40282306e+38  -3.40282306e+38  -3.40282306e+38 ...,  -3.40282306e+38
   -3.40282306e+38  -3.40282306e+38]
 ..., 
 [ -3.40282306e+38  -3.40282306e+38  -3.40282306e+38 ...,  -3.40282306e+38
   -3.40282306e+38  -3.40282306e+38]
 [ -3.40282306e+38  -3.40282306e+38  -3.40282306e+38 ...,  -3.40282306e+38
   -3.40282306e+38  -3.40282306e+38]
 [ -3.40282306e+38  -3.40282306e+38  -3.40282306e+38 ...,  -3.40282306e+38
   -3.40282306e+38  -3.40282306e+38]]
The script took 0.457999944687seconds.
>>> 

我需要做的是随机选择这个数组中的一个位置（单元格），然后在这个点上运行你生成的代码，让洪水填充算法扩展，直到像上面的例子那样达到最大值，或者直到达到设定的单元格数量（用户可以设置洪水填充算法的选择不超过25个单元格）。然后理想情况下，新选中的单元格会作为一个单一的栅格输出，并保持其地理参考结构。

#import modules
from osgeo import gdal
import numpy as np
import os, sys, time

#start timing
startTime = time.time()

#register all of drivers
gdal.AllRegister()

#get raster
geo = gdal.Open("C:/Users/Harmon_work/Desktop/Python_Scratch/all_fill.img")

#read raster as array
arr = geo.ReadAsArray()
data = geo.ReadAsArray(0, 0, geo.RasterXSize, geo.RasterYSize).astype(np.float)
print data

#get image size
rows = geo.RasterYSize
cols = geo.RasterXSize
bands = geo.RasterCount

#get georefrence info
transform = geo.GetGeoTransform()
xOrgin = transform[0]
yOrgin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

#get array dimensions
row = data.shape[0]
col = data.shape[1]

#get random position in array
randx = random.randint(1, row)
randy = random.randint(1, col)
print randx, randy

neighbours = [(-1,-1), (-1,0), (-1,1), (0,1), (1,1), (1,0), (1,-1), (0,-1)]
mask = np.zeros_like(data, dtype = bool)

#start coordinate
stack = [(randx,randy)]

while stack:
    x, y = stack.pop()
    mask[x, y] = True
    for dx, dy in neighbours:
        nx, ny = x + dx, y + dy
        if (0 <= nx < data.shape[0] and 0 <= ny < data.shape[1]
            and not mask[nx, ny] and abs(data[nx, ny] - data[x, y]) <= 1):
            stack.append((nx, ny))

for line in mask:
    print ''.join('01'[i] for i in line)

#run time
endTime = time.time()
print 'The script took ' + str(endTime-startTime) + 'seconds.'

再次感谢你的帮助。如果有什么不清楚的，请问我。

Answer 1

这可以通过一种类似于洪水填充的算法来实现，使用一个栈：

import numpy as np

z = '''33 33 33 37 38 37 43 40
33 33 33 38 38 38 44 40
36 36 36 36 38 39 44 41
35 36 35 35 34 30 40 41
36 36 35 35 34 30 30 41
38 38 35 35 34 30 30 41'''
z = np.array([[int(i) for i in line.split()] for line in z.splitlines()])

neighbours = [(-1,-1), (-1,0), (-1,1), (0,1), (1,1), (1,0), (1,-1), (0,-1)]
mask = np.zeros_like(z, dtype = bool)
stack = [(3,2)] # push start coordinate on stack

while stack:
    x, y = stack.pop()
    mask[x, y] = True
    for dx, dy in neighbours:
        nx, ny = x + dx, y + dy
        if (0 <= nx < z.shape[0] and 0 <= ny < z.shape[1] 
            and not mask[nx, ny] and abs(z[nx, ny] - z[x, y]) <= 1):
            stack.append((nx, ny))

for line in mask:
    print ''.join('01'[i] for i in line)    

结果：

00000000
00000000
11110000
11111000
11111000
00111000