#!/usr/bin/env python
############################################################################
#
# MODULE:    v.nn.py
# AUTHOR(S):	Michael Barton, Arizona State University
# PURPOSE:	Nearest neighbor analysis for 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.
#
# REFERENCES: Kintigh, Keith W. 1990. Intrasite spatial analysis: a 
#             commentary on major methods. In Mathematics and Information 
#             Science in Archaeology: a Flexible Framework, ed. A. Voorrips, 
#             165-200. Bonn: Holos.
#
#             Chen , DongMei and Arthur Getis. 1998. POINT PATTERN 
#             ANALYSIS (PPA). Department of Geography, San Diego State
#             University. http://www.geog.ucsb.edu/~dongmei/ppa/ppa.html 
#             (accessed 1 July 2011)
#
#############################################################################

#%module
#% description: Calculates nearest neighbor coefficients for set of vector points
#% keywords: vector, spatial analysis, statistics, nearest neighbor
#%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

import math
import os

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 NNcalc(map1, map2):
    """calculates nearest neighbor coefficient for a single map
    or for two maps"""
        
    ## Calculate region perimeter and area (rectangular area only)
    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
    perimeter = (2 * ns_extent) + (2 * ew_extent)
    
    ## Calculates 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'])

    ## Calculate expected mean nearest neighbor distance for random distribution of points
    expected_distance = (0.5 * math.sqrt(area/npts2)) + (0.0514 + 0.041/math.sqrt(npts2)) * perimeter / npts2
        
    pointslist = grass.read_command('v.distance', _from=map1, to=map2, 
                                     column='temp', upload='dist', dmin=0.000001, 
                                     flags='p', quiet=True)

    
    pointslist = pointslist.split('\n')
    npoints = len(pointslist)
    
    totaldist = 0.0
    for i in pointslist:
        try:
            val = float(i.split('|')[1])
            totaldist += float(i.split('|')[1])
        except: 
            continue
    
    observed_distance = totaldist / npoints
    return expected_distance / observed_distance

def main():
    """Main method for script"""
    
    map1 = options['input1']
    map2 = options['input2']
    
    
    if map2 == '' or map1 == map2:
        # Single map nearest neighbor for clustering analysis
        map2 = map1
        header = 'Nearest neighbor coeficient for points in %s: ' % map1
    else:
        # Nearest neighbor for correspondence between two point maps
        header = 'Nearest neighbor coeficient of association between points in %s and %s: ' % (map2, map1)
            
    NN = NNcalc(map1, map2)
            
    # Display results in output window
    print '\n%s %G\n' % (header, NN)
    
    
if __name__ == "__main__":
    options, flags = grass.parser()
    main()