python-如何迭代一个numpy数组并选择相邻的单元格

我正在将USGS高程栅格数据集转换为Numpy数组,然后尝试随机选择数组中的位置.从这个位置,我想创建一种方法来识别周围的八个单元,以查看这些单元的高度是否在随机选择的单元的一米以内.

在这里它变得更加复杂…如果邻居在一米之内,则将调用相同的方法,然后重复该过程,直到一米以内不再有像元或所选像元的数量达到规定的限制.

如果不清楚,希望下面的2d数组示例更有意义.随机选择粗体/斜体单元格(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.'

再次感谢你的帮助.如果有任何不清楚的地方,请问我问题.

解决方法:

可以使用类似于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