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

Question

我正在将 USGS 高程栅格数据集转换为 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 脚本界面时遇到了重大问题。

感谢您的帮助！

问题继续……

感谢巴斯斯温克尔的回答。我试图将您的代码合并到我迄今为止编写的代码中，最终得到一个无限循环。下面是我写的。要完成这项工作，我需要克服两个主要步骤。这是从我的栅格生成的数组示例（-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

这可以通过类似于flood fill 的算法使用堆栈来完成：

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