#!/usr/bin/env python
############################################################################
#
# MODULE:    v.lda.py
# AUTHOR(S):	Michael Barton, Arizona State University
# PURPOSE:	Local density analysis of vector points for GRASS 6+
#
# COPYRIGHT:	(C) 2011 by the GRASS Development Team
#
#		This program is free software under the GNU General Public
#		License (>=v2). Read the file COPYING that comes with GRASS
#		for details.
#
# REFERENCE:    Johnson, I. (1984). Cell frequency recording and analysis of 
#    artifact distributions. In Intrasite Spatial Analysis in Archaeology. 
#    H. J. Hietala. Cambridge, Cambridge University Press: 75-96.
#############################################################################

#%module
#% description: Calculates local density coefficients for set of vector points
#% keywords: vector, spatial analysis, statistics, local density analysis
#%End

#%option
#% key: input1
#% type: string
#% gisprompt: old,vector,vector
#% description: Vector points base file
#% required : yes
#%end

#%option
#% key: input2
#% type: string
#% gisprompt: old,vector,vector
#% description: Vector points file to compare with base file
#% required : no
#%end

#%option
#% key: ldafile
#% type: string
#% gisprompt: file,file,file
#% description: File in which to save LDA output
#% required : yes
#% answer: lda.txt
#%end

#%option
#% key: min_r
#% type: double
#% description: Minimum neighborhood radius for local density calculations
#% required : yes
#%end

#%option
#% key: max_r
#% type: double
#% description: Maxium neighborhood radius for local density calculations
#% required : yes
#%end

#%option
#% key: step_r
#% type: double
#% description: Increments to neighborhood radius for local density calculations
#% required : yes
#%end

#%Flag
#%  key: a
#%  description: Append results to existing LDA.txt file instead of overwrite LDA.txt
#%END

#%Flag
#%  key: g
#%  description: Save results in csv format for import into spreadsheet or statistical program
#%END

import math
import os
import subprocess

if not os.environ.has_key("GISBASE"):
    print "You must be in GRASS GIS to run this program."
    sys.exit(1)
    
# GRASS binding
try:
    import grass.script as grass
except ImportError:
    sys.exit(_("No GRASS-python library found"))


def pointdens(radius, ptcount1, ptlist):
    """takes output from v.distance and n of points in map 1 and 
    converts it to mean density of points in
    map 2 around each point of map 1"""
    ptlist = ptlist.split('\n')
    l = []
    neighbors = 0
    for i in ptlist:
        try:
            l.append(int(i.split('|')[0]))
        except:
            continue
    ptcount1 += 1
    for i in range(1, ptcount1):
        neighbors += l.count(i)
        
    # Mean local density of points in map2 around each point in map1 with radius    
    dens = neighbors / (ptcount1 * math.pi * radius**2)
        
    return dens

def main():
    map1 = options['input1']
    map2 = options['input2']
    if map2 == '': map2 = map1
    minradius = float(options['min_r'])
    maxradius = float(options['max_r'])
    step = float(options['step_r'])
    ldafile = options['ldafile']

    ##Gets the number of points in each vector points file
    npts1 = int(grass.parse_command('v.info', map=map1, flags='t')['points'])
    npts2 = int(grass.parse_command('v.info', map=map2, flags='t')['points'])
    totalpts = npts1 * npts2

    ##Calculates region area
    ns_extent = float(grass.parse_command('g.region', flags='eg')['ns_extent'])
    ew_extent = float(grass.parse_command('g.region', flags='eg')['ew_extent'])
    area = ns_extent * ew_extent
    
    # Calculate LDA coefficients for each neighborhood radius
    radius = minradius
    
    #if flags['a'] = True:
#    flags['g'] = False
    
    # Set up text file for LDA output
    if map1 == map2:
        header = 'Local Density Analsis of Clustering of Points in %s' % map1
    else:
        header = 'Local Density Analysis of %s in neighborhood of %s' % (map2, map1)
    
    if flags['a'] == True:
        f = open(ldafile,  'a')
        f.write('\n\n%s' % header)
    else:
        f = open(ldafile, 'w')
        f.write('%s' % header)
        
    if flags['g'] == True:
        f.write('\n"radius","observed","expected","LDA"')
    else:
        f.write('\nradius\tobserved\texpected\tLDA')
        f.write('\n======\t========\t========\t===')
    
    while radius <= maxradius:
        # global density of points in map2
        expected_dens = npts2  / area
        
        # calculate observed density of map2 points around each map1 point
        pointslist = grass.read_command('v.distance', _from=map1, to=map2, 
                                         column='temp', upload='dist', dmax=radius, 
                                         flags='pa', quiet=True)

        observed_dens = pointdens(radius, npts1, pointslist)
        LDA = observed_dens / expected_dens

        if flags['g'] == True:
            f.write('\n%F,%F,%F,%F' % (radius, observed_dens, expected_dens, LDA))
        else:
            f.write('\n%G\t%G\t%G\t%G' % (radius, observed_dens, expected_dens, LDA))
            
        radius += step     

    if flags['g'] == False:
        f.write('\n**********************************************************')
        f.write('\n')
        
    f.close()

    # Display LDA results in output window
    print(" ")
    print("**********************************************************")
    print("Local density coefficients have been written to the file")
    print("%s and are displayed below" %ldafile)
    print("**********************************************************")
    print("")
    f = open(ldafile, 'r')
    print(f.read())

if __name__ == "__main__":
    options, flags = grass.parser()
    main()