#! /usr/bin/env python
# 
""""
Build a kernel for ESRI foacl statistics.
ESRI's focalstatistics can calculate the mean of values
within a given radius of each cell,
but it does not implicitly allow for distance weighting.
However, one can implement wieghting by creating a
kernel file of weights.  It can be tedious to create this
file by hand, so this utility does it for you.
Arguments are optional.  It creates wk[radius].txt.
It then checks the file to make sure it is normalized to 1.0.
See also kernel_annulus.py at
http://gis.ess.washington.edu/projects/zcolor
Harvey Greenberg hgreen@uw.edu Mon, Oct 31, 2016 12:31:49 PM
"""""
print "kernel.py verserion 1.0"
import numpy as np
import sys
if len(sys.argv) < 2:
  print "No first argument for radius.  Using default."
  radius = 4
else:
  radius = int(sys.argv[1])
if len(sys.argv) < 3:
  print "No second argument for exponent.  Using default."
  exponent = 1    # weight proportional to distance
else:
  exponent = float(sys.argv[2])
print "Radius is",radius,"cells, exponent is",exponent
exponent /= 2.  # build in the square root for the hypotenuse.
edge = 2*radius + 1
box = np.array(edge*edge*[0.], dtype='f4').reshape(edge,edge)
for i in range(radius+1):
  for j in range(radius+1):
    if (i==0) and (j==0):
      weight = 2.
    else:
      weight = 1. / ((i**2 + j**2)**exponent)
    box[radius+i,radius+j] = weight
    box[radius+i,radius-j] = weight
    box[radius-i,radius+j] = weight
    box[radius-i,radius-j] = weight

total = box.sum() / (edge*edge)	# fixed Jan 11
for i in range(edge):
  for j in range(edge):
    box[i,j] /= total


f=open("wk"+`radius`+".txt","w")
f.write("%i %i\n" % (edge,edge))
for i in range(edge):
  for j in range(edge):
    f.write("%f " % box[i,j])
  f.write("\n")

f.close()
# check rounding error:
summ = 0.
f=open("wk"+`radius`+".txt")
f.readline()
for i in range(edge):
  nums = f.readline().split()
  for val in nums:
    summ += float(val)

print summ / (edge*edge), "on wk"+`radius`+".txt"