#-------------------------------------------------------------------------------
# Name:			Sediment Wedge Model Generator Aug 2013 Part V.py
# Purpose:
#
#
# Author:		Dr. Arghya Goswami
# Contact:		arghya.goswami@gmail.com
#
# Created:		06/08/2013
#
# Copyright:	(c) Dr. Arghya Goswami 2013
#
# Licence:		This program is free software; you can redistribute it and/or
#				modify it under the terms of the GNU General Public License,
#				version 2, as published by the Free Software Foundation.
#
#				This program is distributed in the hope that it will be useful,
#				but WITHOUT ANY WARRANTY; without even the implied warranty of
#				MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#				GNU General Public License for more details.
#
#				You should have received a copy of the GNU General Public License
#				along with this program; if not, write to:
#
#				Free Software Foundation, Inc.,
#				51 Franklin St, Fifth Floor,
#				Boston, MA, 02110-1301  USA
#-------------------------------------------------------------------------------
#!/usr/bin/env python
#
#-------------------------------------------------------------------------------
import arcpy, sets, numpy, sys, os, shutil, glob, time, fileinput, string, Tkinter, tkFileDialog
from operator import itemgetter
from Tkinter import *
from arcpy.sa import *
"""-----------------------------------------------------------------------------
# Start Time
-----------------------------------------------------------------------------"""
StartTime = time.asctime(time.localtime(time.time()))
print "Start Time: ", StartTime,"\n"
"""-----------------------------------------------------------------------------
# Check Script name
-----------------------------------------------------------------------------"""
#scriptname = os.path.basename(__file__)
#if "module" in scriptname:
#	print "Please SAVE the script with a Meaningful NAME and RE-RUN!\n"
"""-----------------------------------------------------------------------------
# Definitions used in this Script
-----------------------------------------------------------------------------"""
""" Browse a file function """
def BrowseFile(msg1,indir,ext):
	#Creating Main Window
	root = Tk()
	if ext == "csv":
		desc1 = "Comma Seperated Value Files"
		desc2 = "*." + ext
	elif ext == "txt":
		desc1 = "Text Files"
		desc2 = "*." + ext
	elif ext == "shp":
		desc1 = "Shape Files"
		desc2 = "*." + ext
    # Creating Secondary Window for File Browsing
	filename = tkFileDialog.askopenfilename(filetypes = [(desc1,desc2),('All files','*')], title = msg1, initialdir = indir)
    # Close the Main Window
	root.withdraw()
    # Return Filename
	return filename

""" Browse a directory function """
def BrowseDir(msg2,indir):
    #Creating Main Window
	root = Tk()
    # Creating Secondary Window for Folder Browsing
	dirname = tkFileDialog.askdirectory(title = msg2, initialdir=indir)
    # Close the Main Window
	root.withdraw()
    # Return Dirname
	return dirname

""" Check the existance of a directory function """
def DirCheck(filename):
	# Check if the string is a file name or directory name
	if filename[len(filename)-4] == ".":
		dname = os.path.dirname(filename)
	else:
		dname = filename
	# Check if Dir exists, else create
	if os.path.exists(dname):
	    pass
	else:
	    os.mkdir(dname)
	# Check if the path name if Dir ends with a '/' else add it
	if dname[(len(dname)-1):]=="\\" or dname[(len(dname)-1):] == "/":
		dpath = dname
	else:
		dpath = dname + "/"
	return dpath

""" Returns field names in a shapefile """
def AllFields(infile):
	FieldNnames = []
	fields = arcpy.ListFields(infile)
	FieldNnames = []
	for field in fields:
		FieldNnames.append(field.name)
	return FieldNnames

""" Create the list of files in a folder """
def getFileList(root,ext):
	# for folder path either use "Browse"/"browse" or provide the path to a folder
	# for all filetypes send extention as an empty string (""), else use like ".csv"
	filelist = []
	if root == "Browse" or root == "browse":
		root = BrowseDir("Please select the Directory", "C:/")
	# select the type of file, for instance *.jpg or all files *.*
	for folder in glob.glob(root):
		print "Files in", folder, "is being populated to a list"
	if ext == "":
		ext = "/*.*"
	else:
		ext = "/*" + ext
	for file in glob.glob(folder + ext):
	# retrieves the stats for the current file as a tuple
	# (mode, ino, dev, nlink, uid, gid, size, atime, mtime, ctime)
	# the tuple element mtime at index 8 is the last-modified-date
		filelist.append(file)
	return filelist
	print "\nFile List Prepared, starting to list them......\n"
"""-----------------------------------------------------------------------------
# Main Script
-----------------------------------------------------------------------------"""
""" Step 11: Calculate Averages distance for D2WS & D2EN Fields """

sid = [17] #[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18]

for sids in sid:
	# Select Input Merged PseudoPoint file
	inFC = "C:/OES/PaleoBathymetry/Projected Shapes/diced/" + ((str(sids)).zfill(2) + "/Done/" + (str(sids)).zfill(2) + "_Merged_PseudoPoints_V3_RajBathy.V3.shp").replace(" ","") # BrowseFile("Select Merged Shapefile","C:/OES/PaleoBathymetry/Projected Shapes/diced/","shp")
	print "Processing ",inFC

	# Export to CSV
	fieldNames = ["SetID","LineID","WSTx","WSTy","WSCName","ENTx","ENTy","ENCName","LON","LAT","ET1","ET5","D2WS","D2EN","DST000","MMODSED000","RMODBATHY0","MMODB000"]
	outFCcsv = inFC.replace(".shp",".csv")
##	arcpy.ExportXYv_stats (inFC, fieldNames, "COMMA", outFCcsv, "ADD_FIELD_NAMES")

	# Make Shapefile with original lat & lon
##	try:
##		t = t+1
##	except:
##		t =0
##	out_Layer = "XYEventLayer" + str(t)
##	spref = "Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj"
##	arcpy.MakeXYEventLayer_management(outFCcsv, "LON", "LAT", out_Layer, spref)
	outFCshp2 = inFC.replace("shp","shp")
##	outFCshp2_folder = os.path.dirname(outFCshp2)
##	outFCshp2_shpfile = os.path.basename(outFCshp2)
##	arcpy.FeatureClassToFeatureClass_conversion(out_Layer, outFCshp2_folder, outFCshp2_shpfile)
##	arcpy.DeleteField_management(outFCshp2,"XCoord")
##	arcpy.DeleteField_management(outFCshp2,"YCoord")
	print "V2 Created....\n"

##	# Spatially Join the V2 PseudoPoint to Divine2003 data
##	join_features = "C:/OES/PaleoBathymetry/Divine2003/Divine_2003_ST.shp"
##	outFCshp3 = inFC.replace("shp","V3.shp")
##	arcpy.SpatialJoin_analysis (outFCshp2, join_features, outFCshp3, "JOIN_ONE_TO_ONE", "KEEP_COMMON", "", "CLOSEST","","Distance")
##	arcpy.DeleteField_management(outFCshp3,"Join_Count")
##	arcpy.DeleteField_management(outFCshp3,"TARGET_FID")
##	print "V3 Created....\n"

	# Get the distance fields and make output file
	field1 = "D2WS" # raw_input("Enter First Distance Field Name")
	rl = "Selected first distance field is " + field1 + "\n"
	##fReport.write(rl)
	print rl

	field2 = "D2EN" # raw_input("Enter Second Distance Field Name")
	rl = "Selected second distance field is " + field2 + "\n"
	##fReport.write(rl)
	print rl

	outF = DirCheck(outFCshp2) + field1 + "-" + field2 + "_Distace_Averaged_withRajBathySed1.csv"
	rl = "Output file is " + outF + "\n"
	##fReport.write(rl)
	print rl

	fout= open(outF,'w')
##	header = "Distance,\n"
	header = "Distance,WSAvgST,WSAvgET5,WSAvgET1,WSAvgRB,WSAvgMB,WSAvgMS,WSAvgRST,ENAvgST,ENAvgET5,ENAvgET1,ENAvgRB,ENAvgMB,ENAvgMS,ENAvgRST\n"
	fout.write(header)

	d = 0.0
	offset = 1000
	# Loop through increasing distance until distace = 2000 Km
	while d <=2500000:  # 3000 Km in meters
##		# Show Progress
##		if (int(d)%33000==0):
##			print str(d/25000), " Km processed \n"
##		print d

		# Reset variables
		WSRSTvals = []    # list with ST vals from West End
		WSSTvals = []    # list with ST vals from West End
		WSET5vals = []    # list with ET5 vals from West End
		WSET1vals = []    # list with ET1 vals from West End
		WSRBvals = []
		WSMBvals = []
		WSMSvals = []

		ENRSTvals = []    # list with ST vals from East End
		ENSTvals = []    # list with ST vals from East End
		ENET5vals = []    # list with ET5 vals from East End
		ENET1vals = []    # list with ET5 vals from East End
		ENRBvals = []
		ENMBvals = []
		ENMSvals = []

		WSAvgRST = 0.0
		WSAvgST = 0.0
		WSAvgET5 = 0.0
		WDAvgET1 = 0.0
		WSAvgRB = 0.0
		WSAvgMB = 0.0
		WSAvgMS = 0.0

		ENAvgRST = 0.0
		ENAvgST = 0.0
		ENAvgET5 = 0.0
		ENAvgET1 = 0.0
		ENAvgRB = 0.0
		ENAvgMB = 0.0
		ENAvgMS = 0.0

		# Where_clause for "DfWKm" and related calculations
		where_clause_ws = "\"" + field1 + "\" >= "+ str(d-1000) + "AND \"" + field1 + "\" <= " + str(d+1000)
		# print where_clause_w

		WSRSTvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "DST000", where_clause_ws)]
		# print len(WSTvals)
		# print WSTvals
	##	del row
		if len(WSRSTvals)>1:
			WSAvgRST = numpy.average(WSRSTvals)
		elif len(WSRSTvals)==0:
			WSAvgRST = ""
		else:
			 WSAvgRST = WSRSTvals[0]
		# print WSAvgRST

		WSSTvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "DST000", where_clause_ws)]
		# print len(WSTvals)
		# print WSTvals
	##	del row
		if len(WSSTvals)>1:
			WSAvgST = numpy.average(WSSTvals)
		elif len(WSSTvals)=="":
			WSAvgST = 0
		else:
			 WSAvgST = WSSTvals[0]
		# print WAvgST
		WSET5vals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "ET5", where_clause_ws)]
		# print len(WET5vals)
		# print WET5vals
	##	del row
		if len(WSET5vals)>1:
			WSAvgET5 = numpy.average(WSET5vals)
		if len(WSET5vals)== "":
			WSAvgET5 = 0
		else:
			WSAvgET5 = WSET5vals[0]
		# print WAvgET5
		WSET1vals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "ET1", where_clause_ws)]
		# print len(WSET1vals)
		# print WSET1vals
	##	del row
		if len(WSET1vals)>1:
			WSAvgET1 = numpy.average(WSET1vals)
		if len(WSET1vals)=="":
			WSAvgET1 = 0
		else:
			WSAvgET1 = WSET1vals[0]
		# print WAvgET1

		WSRBvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "RMODBATHY0", where_clause_ws)]
		# print len(WSRBvals)
		# print WSRBvals
	##	del row
		if len(WSRBvals)>1:
			WSAvgRB = numpy.average(WSRBvals)
		if len(WSRBvals)=="":
			WSAvgRB = 0
		else:
			WSAvgRB = WSRBvals[0]
		# print WSAvgRB

		WSMBvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "MMODB000", where_clause_ws)]
		# print len(WSMBvals)
		# print WSMBvals
	##	del row
		if len(WSMBvals)>1:
			WSAvgMB = numpy.average(WSMBvals)
		if len(WSMBvals)=="":
			WSAvgMB = 0
		else:
			WSAvgMB = WSMBvals[0]
		# print WSAvgMB

		WSMSvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "MMODSED000", where_clause_ws)]
		# print len(WSMSvals)
		# print WSMSvals
	##	del row
		if len(WSMSvals)>1:
			WSAvgMS = numpy.average(WSMSvals)
		if len(WSMSvals)=="":
			WSAvgMS = 0
		else:
			WSAvgMS = WSMSvals[0]
		# print WSAvgMS


		# Where_clause for "DfEKm" and related calculations
		where_clause_en = "\"" + field2 + "\" >= "+ str(d-1000) + "AND \"" + field2 + "\" <= " + str(d+1000)
		# print where_clause_en

		ENRSTvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "DST000", where_clause_en)]
		# print len(ENRSTvals)
		# print ENRSTvals
	##	del row
		if len(ENRSTvals)>1:
			ENAvgRST = numpy.average(ENRSTvals)
		if len(ENRSTvals)=="":
			ENAvgRST = 0
		else:
			ENAvgRST = ENRSTvals[0]
		# print ENAvgRST

		ENSTvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "DST000", where_clause_en)]
		# print len(ESTvals)
		# print ESTvals
	##	del row
		if len(ENSTvals)>1:
			ENAvgST = numpy.average(ENSTvals)
		if len(ENSTvals)=="":
			ENAvgST = 0
		else:
			ENAvgST = ENSTvals[0]
		# print EAvgST
		ENET5vals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "ET5", where_clause_en)]
		# print len(EET5vals)
		# print EET5vals
	##	del row
		if len(ENET5vals)>1:
			ENAvgET5 = numpy.average(ENET5vals)
		if len(ENET5vals)== "":
			ENAvgET5 = 0
		else:
			ENAvgET5 = ENET5vals[0]
		# print EAvgET5
		ENET1vals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "ET1", where_clause_en)]
		# print len(EET1vals)
		# print EET1vals
	##	del row
		if len(ENET1vals)>1:
			ENAvgET1 = numpy.average(ENET1vals)
		if len(ENET1vals)== "":
			ENAvgET1 = 0
		else:
			ENAvgET1 = ENET1vals[0]
		# print EAvgET1


		ENRBvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "RMODBATHY0", where_clause_en)]
		# print len(ENRBvals)
		# print ENRBvals
	##	del row
		if len(ENRBvals)>1:
			ENAvgRB = numpy.average(ENRBvals)
		if len(ENRBvals)=="":
			ENAvgRB = 0
		else:
			ENAvgRB = ENRBvals[0]
		# print ENAvgRB

		ENMBvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "MMODB000", where_clause_en)]
		# print len(ENMBvals)
		# print ENMBvals
	##	del row
		if len(ENMBvals)>1:
			ENAvgMB = numpy.average(ENMBvals)
		if len(ENMBvals)=="":
			ENAvgMB = 0
		else:
			ENAvgMB = ENMBvals[0]
		# print ENAvgMB

		ENMSvals = [row[0] for row in arcpy.da.SearchCursor(outFCshp2, "MMODSED000", where_clause_en)]
		# print len(ENMSvals)
		# print ENMSvals
	##	del row
		if len(ENMSvals)>1:
			ENAvgMS = numpy.average(ENMSvals)
		if len(ENMSvals)== "":
			ENAvgMS = 0
		else:
			ENAvgMS = ENMSvals[0]
		# print ENAvgMS


		# Create line and add to output file
##		wLine = str(d)+ "\n"
				#Distance,			WSAvgST,			WSAvgET5,				WSAvgET1,			WSAvgRB,			WSAvgMB,			WSAvgMS,				ENAvgST,			ENAvgET5,				ENAvgET1,			ENAvgRB,			ENAvgMB,			ENAvgMS
		wLine = str(d) + "," + str(WSAvgST) + "," + str(WSAvgET5)  + "," + str(WSAvgET1) + "," + str(WSAvgRB) + "," + str(WSAvgMB) + "," + str(WSAvgMS) + "," + str(WSAvgRST) + "," + str(ENAvgST) + "," + str(ENAvgET5) + "," + str(ENAvgET1) + ","+ str(ENAvgRB) + "," + str(ENAvgMB) + ","+ str(ENAvgMS) + "," + str(ENAvgRST) + "\n"
		fout.write(wLine)

		# increment in distance (d) value
		d = d + offset
"""-----------------------------------------------------------------------------
# End Time
-----------------------------------------------------------------------------"""
EndTime = time.asctime(time.localtime(time.time()))
print "\nEnd Time: ", EndTime,"\n"
