from scipy.spatial import cKDTree
import numpy as np
from dbfpy import dbf
import cmath

def gaussW(dists,band):
    out = np.zeros(dists.shape)
    for i in xrange(0,len(out)):
        temp = np.power(dists[i],2)/(2.0*np.power(float(band),2))
        out[i] = np.exp(-1.0 * temp)          
    return out

filestring1 = #Grid data filestring
filestring2 = #Candidate Data Point filestring

pat = dbf.Dbf(filestring2)
grid = dbf.Dbf(filestring1)

treepoints = np.zeros((len(pat),2))
angles = np.zeros((len(pat)))

for i in xrange(0,len(pat)):
    rec = pat[i]
    treepoints[i][0] = rec['POINT_X']
    treepoints[i][1] = rec['POINT_Y']
    angles[i] = rec['DirectUsed']
    if i % 10000 == 0:
        print i

tree = cKDTree(treepoints)

# in tree, k is artificially high (300000) and bandwidth is final value (100)
output = []
for i in xrange(0,len(grid)):
    rec = grid[i]
    testpoint = [[rec['X_Centroid'],rec['Y_Centroid']]]
    res,idx = tree.query(testpoint,300000,0,2,100)
    res = res[0][np.where(res[0] < np.Inf)[0]]
    if len(res) > 0:
        idx = idx[0][:len(res)]
        weight = gaussW(res,100)
        thetas = angles[idx]
        # create list of complex numbers
        cThetas = []
        for i in xrange(0,len(thetas)):
            cThetas.append(complex(np.cos(thetas[i]),np.sin(thetas[i])))
        cThetas = np.array(cThetas)
        #calculate weighted mean direction
        temp = np.sum(weight*cThetas)/np.absolute(np.sum(weight*cThetas))
        rec['MeanDirect'] = np.angle(temp,deg=True)
        rec.store()
    else:
        # If there are no local points, record -999 to exclude later.
        rec['MeanDirect'] = -999
        rec.store()
    if i % 1000 == 0:
        print i
pat.close()
grid.close()