Using Numpy to generate ASCII kernel

Question

I am using ArcGIS focal statistics tool to add spatial autocorrelation to a random raster to model error in DEMs. The input DEM has a 1.5m pixel size and the semivariogram exhibits a sill around 2000m. I want to make sure to model the extent of the autocorrelation in the input in my model.

Unfortunately, ArcGIS requires that the input kernel be in ASCII format, where the first line defines the size and the subsequent lines define the weights.

Example:

I need to generate a 1333x1333 kernel with an inverse distance weighting and immediately went to python to get this done. Is it possible to generate a matrix in numpy and assign values by ring? Does a better programmatic tool within numpy exist to generate a plain text matrix.

This is similar to this question, but I need to have a fixed central value and descending rings, as per the example above.

Note: I am a student, but this not a homework assignment...those ended years ago. This is a part of a larger research project that I am working on and any help (even just a nudge in the right direction) would be appreciated. The focus of this work is not programming kernels, but exploring errors in DEMs.

Just as a warning, I'd be suprised if Arc can handle focal statistics with a 1333x1333 kernel... Good luck, at any rate! @DSM & @wim already have your answer, but if you run into further problems, feel free to ask. You'll find that numpy (and particularly scipy.ndimage with what you're doing) can be very capable if you understand how to "glue" a few lower-level functions together to do what you want. That having been said, the python community is still lacking a good kriging/flexible-interpolation-of-any-sort library, at the moment. — Joe Kington, Commented Mar 5, 2012 at 2:32
@JoeKington Going to give the gigantic kernel a try to see what it can do. If Arc fails, I might give it a go through numpy. If I can use numpy to generate a spatially autocorrelated DEM with a standard deviation of 1 I will be all set to go. Will report back! — Jzl5325, Commented Mar 5, 2012 at 11:50
Currently running the focal statistics portion of the model...49037% done on a 6MB ASCII kernel...looks like I will be implementing this step using numpy. — Jzl5325, Commented Mar 6, 2012 at 2:01

wim · Accepted Answer · 2012-03-05 02:27:23Z

3

I'm not sure if there is a built-in way, but it should not be hard to roll your own:

>>> 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

edited Mar 5, 2012 at 2:27

answered Mar 5, 2012 at 1:46

wim

365k110 gold badges669 silver badges798 bronze badges

Thanks! To make sure I understand what is going on: In the first def statement you are creating an empty array of zeros which is the correct size. Then you iterate through (xrange) using numpy's slicing(?) adding one with each iteration and appending the array of zeros? The second def sets the shape, writes the first line, and converts to string? Again thanks!
– Jzl5325
Commented Mar 5, 2012 at 11:50
pretty much.. you can think of the loop as working like this: the kernel is made up of a sum of stack of "squares" of ones, and the squares are decreasing in area - each time in the loop is adding another square onto the middle of the stack.
– wim
Commented Mar 5, 2012 at 13:34
Thanks. To add a note for people running python 2.6 or less. In the kernel_to_string return statement be sure to include positional arguments, ie {0} {1}.
– Jzl5325
Commented Mar 6, 2012 at 1:42

Add a comment |

DSM · Accepted Answer · 2012-03-05 02:27:01Z

[Hmmph. @wim beat me, but I'd already written the following, so I'll post it anyway.] Short version:

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")

which produces

~/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 explanation of the magic: the first thing we're going to need are the indices of each entry in the matrix, and for that we can use 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]])

So, taken pairwise, these are the x and y coordinates of each element of a 5x5 array. The centre will be at N//2,N//2 (where // is truncating division), so we can subtract that to get the distances, and take the absolute value because we don't care about the sign:

>>> 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]])

Now looking at the original grid, it looks like you want the maximum value of the two:

>>> 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]])

which looks good but goes the wrong way. We can invert, though:

>>> 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]])

and you want 1-indexing, it looks like, so:

>>> 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]])

and there we have it. Everything else is boring IO.

pretty cool.. i was looking for a tricky way to do it with np.tri, np.diag, maybe even a convolution, before i settled on the simple broadcasting loop :) easy to get distracted by "interesting" solutions on things like these.. — wim, Commented Mar 5, 2012 at 2:33

Collectives™ on Stack Overflow

Using Numpy to generate ASCII kernel

2 Answers 2

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Linked

Related