#! /usr/bin/env python
"""" 
Build a kernel for ESRI foacl statistics.
ESRI's focalstatistics can calculate the mean of values
within a given radius (or within an annulus) 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 annulus[radius].txt.
It then checks the file to make sure it is normalized to 1.0.
See also kernel.py at
http://gis.ess.washington.edu/projects/zcolor
Harvey Greenberg hgreen@uw.edu Mon, Oct 31, 2016 12:31:49 PM
""""" 
print "kernel_annulus.py version 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 inside radius.  Using default."
  inside = radius * .4
else:
  inside = radius * float(sys.argv[2])

if len(sys.argv) < 4:
  print "No third argument for outside radius.  Using default."
  outside = radius * .7
else:
  outside = radius * float(sys.argv[3])

if len(sys.argv) < 5:
  print "No fourth argument for exponent.  Using default."
  exponent = 1    # weight proportional to distance
else:
  exponent = float(sys.argv[4])
print "Total radius is",radius,"cells, Inside radius is",inside,"cells, Outside radius is",outside,"exponent is",exponent
### exponent /= 2.  # build in the square root for the hypotenuse.
edge = 2*radius +1
#rad = 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):
    dist = ((i**2 + j**2)**.5)
    if dist < inside:
      weight = inside-dist
    elif dist > outside:
      weight = dist-outside
    else:
      weight = 1.
    if weight < 1.:
      weight = 1.
    print i,j,dist,weight   
    weight = 1. / (weight**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("annulus"+`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("annulus"+`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 annulus"+`radius`+".txt"	# 1.00003