Использование Numpy для генерации ядра ASCII

Question

Я использую инструмент фокусной статистики ArcGIS для добавления пространственной автокорреляции к случайному растру для моделирования ошибок в ЦМР.Входная ЦМР имеет размер пикселя 1,5 м, а вариограмма имеет порог около 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, чтобы это сделать.Можно ли сгенерировать матрицу в numy и присвоить значения по кольцу?Существует ли лучший программный инструмент в numpy для генерации матрицы простого текста.

Это похоже на этот вопрос , но мне нужно иметь фиксированное центральное значение и нисходящие кольца, как указано впример выше.

Примечание: я студент, но это не домашнее задание ... оно закончилось много лет назад.Это часть более крупного исследовательского проекта, над которым я работаю, и любая помощь (даже просто толчок в правильном направлении) будет принята с благодарностью.Основное внимание в этой работе уделяется не программированию ядер, а изучению ошибок в ЦМР.

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

---

Verboseобъяснение магии: первое, что нам понадобится, - это индексы каждой записи в матрице, и для этого мы можем использовать 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]])

Итак, взялипопарно, это координаты x и y каждого элемента массива 5x5.Центр будет в 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-индексирование, это выглядит так:

>>> 1 + N//2 - numpy.maximum(abs(xx-N//2), abs(yy-N//2))
array([[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]])

и вот оно у нас.Все остальное скучно в IO.