使用Numpy生成ASCII内核

Question

我正在使用ArcGIS焦点统计工具将空间自相关添加到随机栅格来模拟DEM中的错误。输入的DEM具有150万像素的尺寸，半变异函数在2000米左右呈现出一个基石。我想确保在模型中对输入中的自相关范围进行建模。

不幸的是，ArcGIS要求输入内核采用ASCII格式，其中第一行定义了大小，后面的行定义了权重。

例子：

5 5 
1 1 1 1 1 
1 2 2 2 1 
1 2 3 2 1 
1 2 2 2 1 
1 1 1 1 1 

我需要生成与距离倒数加权一个1333x1333的内核，并立即前往Python来完成这件事。是否有可能在numpy中生成矩阵并通过ring来分配值？在numpy中是否存在一个更好的编程工具来生成纯文本矩阵。

这与this question类似，但我需要有一个固定的中心值和降环，如上例所示。

注：我是一名学生，但这不是一项家庭作业......那些年前才结束。这是我正在研究的一个更大的研究项目的一部分，任何帮助（甚至只是在正确的方向推动）都将不胜感激。这项工作的重点不是编程内核，而是探索DEM中的错误。

Answer 1

我不知道是否有一个内置的方式，但它不应该是很难推出自己的：

>>> def kernel_thing(N): 
... import numpy as np 
... n = N // 2 + 1 
... a = np.zeros((N, N), dtype=int) 
... for i in xrange(n): 
...  a[i:N-i, i:N-i] += 1 
... return a 
... 
>>> def kernel_to_string(a): 
... return '{} {}\n'.format(a.shape[0], a.shape[1]) + '\n'.join(' '.join(str(element) for element in row) for row in a) 
... 
>>> print kernel_to_string(kernel_thing(5)) 
5 5 
1 1 1 1 1 
1 2 2 2 1 
1 2 3 2 1 
1 2 2 2 1 
1 1 1 1 1 
>>> print kernel_to_string(kernel_thing(6)) 
6 6 
1 1 1 1 1 1 
1 2 2 2 2 1 
1 2 3 3 2 1 
1 2 3 3 2 1 
1 2 2 2 2 1 
1 1 1 1 1 1 
>>> print kernel_to_string(kernel_thing(17)) 
17 17 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 
1 2 3 3 3 3 3 3 3 3 3 3 3 3 3 2 1 
1 2 3 4 4 4 4 4 4 4 4 4 4 4 3 2 1 
1 2 3 4 5 5 5 5 5 5 5 5 5 4 3 2 1 
1 2 3 4 5 6 6 6 6 6 6 6 5 4 3 2 1 
1 2 3 4 5 6 7 7 7 7 7 6 5 4 3 2 1 
1 2 3 4 5 6 7 8 8 8 7 6 5 4 3 2 1 
1 2 3 4 5 6 7 8 9 8 7 6 5 4 3 2 1 
1 2 3 4 5 6 7 8 8 8 7 6 5 4 3 2 1 
1 2 3 4 5 6 7 7 7 7 7 6 5 4 3 2 1 
1 2 3 4 5 6 6 6 6 6 6 6 5 4 3 2 1 
1 2 3 4 5 5 5 5 5 5 5 5 5 4 3 2 1 
1 2 3 4 4 4 4 4 4 4 4 4 4 4 3 2 1 
1 2 3 3 3 3 3 3 3 3 3 3 3 3 3 2 1 
1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

Answer 2

[Hmmph。 @wim打我，但我已经写了下面的，所以无论如何，我会张贴]短版：

import numpy 

N = 5 

# get grid coords 
xx, yy = numpy.mgrid[0:N,0:N] 
# get the distance weights 
kernel = 1 + N//2 - numpy.maximum(abs(xx-N//2), abs(yy-N//2)) 

with open('kernel.out','w') as fp: 
    # header 
    fp.write("{} {}\n".format(N, N)) 
    # integer matrix output 
    numpy.savetxt(fp, kernel, fmt="%d") 

产生的神奇

~/coding$ python kernel.py 
~/coding$ cat kernel.out 
5 5 
1 1 1 1 1 
1 2 2 2 1 
1 2 3 2 1 
1 2 2 2 1 
1 1 1 1 1 

---

冗长的解释：我们将首先需要的是矩阵中的每个条目的指标，为此，我们可以使用mgrid：

>>> import numpy 
>>> N = 5 
>>> xx, yy = numpy.mgrid[0:N,0:N] 
>>> xx 
array([[0, 0, 0, 0, 0], 
     [1, 1, 1, 1, 1], 
     [2, 2, 2, 2, 2], 
     [3, 3, 3, 3, 3], 
     [4, 4, 4, 4, 4]]) 
>>> yy 
array([[0, 1, 2, 3, 4], 
     [0, 1, 2, 3, 4], 
     [0, 1, 2, 3, 4], 
     [0, 1, 2, 3, 4], 
     [0, 1, 2, 3, 4]]) 

因此，采取配对方式，这些是5x5阵列的每个元素的x和y坐标。中心位于N // 2，N // 2（其中//是截断除法），因此我们可以减去该距离以获得距离，并取绝对值，因为我们不关心该符号：

>>> abs(xx-N//2) 
array([[2, 2, 2, 2, 2], 
     [1, 1, 1, 1, 1], 
     [0, 0, 0, 0, 0], 
     [1, 1, 1, 1, 1], 
     [2, 2, 2, 2, 2]]) 
>>> abs(yy-N//2) 
array([[2, 1, 0, 1, 2], 
     [2, 1, 0, 1, 2], 
     [2, 1, 0, 1, 2], 
     [2, 1, 0, 1, 2], 
     [2, 1, 0, 1, 2]]) 

现在看看原来的网格，它看起来像你想两者的最大值：

>>> numpy.maximum(abs(xx-N//2), abs(yy-N//2)) 
array([[2, 2, 2, 2, 2], 
     [2, 1, 1, 1, 2], 
     [2, 1, 0, 1, 2], 
     [2, 1, 1, 1, 2], 
     [2, 2, 2, 2, 2]]) 

这看起来不错，但去了错误的方式。我们可以反转，虽然：

>>> N//2 - numpy.maximum(abs(xx-N//2), abs(yy-N//2)) 
array([[0, 0, 0, 0, 0], 
     [0, 1, 1, 1, 0], 
     [0, 1, 2, 1, 0], 
     [0, 1, 1, 1, 0], 
     [0, 0, 0, 0, 0]]) 

，你想1索引，它看起来像，所以：

​​

和我们有它。其他一切都是无聊的IO。