#-------------------------------------------------------------------------------
# Name:			SWMGenerator Part II.py
# Purpose:      Sediment Wedge Model Generator
#
#
# 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 AllfFields(infile):
	global 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"
	global 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)
	print "\nFile List Prepared, starting to list them......\n"
"""-----------------------------------------------------------------------------
# Main Script
-----------------------------------------------------------------------------"""
# Global variable
filelist = []
fieldNames = []

# Local Variables

""" Generate a Report file for the script """
# Note: for every run of the script a report file will be generated in "C:/Temp/"
report = DirCheck("C:/Temp/") + (str(StartTime)).replace(":","-") + ".txt"
fReport = open(report,'w')
rl = "Start Time: " + str(StartTime) + "\n"
fReport.write(rl)

"""Step: 6:
    A) Add start lat (sLat), start lon (sLon), end lat (eLat) & end lon (eLon), plus delete "ID" for all Normal to XSLine shapefile in a directory
    B) Join the Nomal2XSLine shapefile to corrosponding ETpoint file"""
root = BrowseDir("Select the forlder with Nomal2XSLine shapefiles: ","C:/OES/")
filelist = []
# 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.\n"
	#raw_input()

for file in glob.glob(folder + '/*Normals2XSectionLine.shp'):

	# 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)

print "\nFile List Prepared, starting conversion......"

for items in filelist:
	items = items.replace("\\","/")
	print "Processing ", items, "..."
	rl = "Total number of features in " + items + " = " + str(arcpy.GetCount_management (items)) + "\n"
	fReport.write(rl)
	print rl
	del rl
	# Add sLon, sLat, eLon, eLat & delete ID
	arcpy.AddField_management (items, "sLon", "DOUBLE")
	arcpy.CalculateField_management(items, "sLon", "!Shape!.firstPoint.X","PYTHON_9.3","")
	arcpy.AddField_management (items, "sLat", "DOUBLE")
	arcpy.CalculateField_management(items, "sLat", "!Shape!.firstPoint.Y","PYTHON_9.3","")
	arcpy.AddField_management (items, "eLon", "DOUBLE")
	arcpy.CalculateField_management(items, "eLon", "!Shape!.lastPoint.X","PYTHON_9.3","")
	arcpy.AddField_management (items, "eLat", "DOUBLE")
	arcpy.CalculateField_management(items, "eLat", "!Shape!.lastPoint.Y","PYTHON_9.3","")
	arcpy.DeleteField_management (items, "ID")
	print "Start - end lon and lat added and 'ID' field deleted in ", items

	# Spatial join of ETPoints to Normals2XSectionLine
	outF6B = items.replace(".shp","_wET.shp")
	arcpy.SpatialJoin_analysis ((items.replace("Normals2XSectionLine","LineETPoints")),items, outF6B)
	print "ETPoints to Normals2XSectionLine join complete...."
	print "Deleting unnecessary fields...."
	arcpy.DeleteField_management (outF6B, "Distance")
	print "'Distance' deleted..."
	arcpy.DeleteField_management (outF6B, "Join_Count")
	print "'Join_Count' deleted..."
	arcpy.DeleteField_management (outF6B, "TARGET_FID")
	print "'TARGET_FID' deleted..."
	arcpy.AddField_management (outF6B, "pLon", "DOUBLE")
	arcpy.CalculateField_management(outF6B, "pLon", "!sLon!","PYTHON_9.3","")
	arcpy.AddField_management (outF6B, "pLat", "DOUBLE")
	arcpy.CalculateField_management(outF6B, "pLat", "!sLat!","PYTHON_9.3","")
	arcpy.DeleteField_management (outF6B, "sLon")
	print "'sLon' deleted..."
	arcpy.DeleteField_management (outF6B, "sLat")
	print "'sLat' deleted..."

	# extract line info from XSectionLine shapefile
	outF6C = items.replace("Normals2XSectionLine","XSectionLine")
	SC = arcpy.SearchCursor(outF6C)
	for row in SC:
		setid = row.getValue("SetID")
		lineid = row.getValue("LineID")
		wstx = row.getValue("WSTx")
		wsty = row.getValue("WSTy")
		wscn = "\"" + str(row.getValue("WSCName")) + "\""
		entx = row.getValue("ENTx")
		enty = row.getValue("ENTy")
		encn = "\"" + str(row.getValue("ENCName"))+ "\""
		row = SC.next()
	del row,SC
	arcpy.AddField_management (outF6B, "SetID", "SHORT")
	arcpy.CalculateField_management(outF6B, "SetID", setid,"PYTHON_9.3","")

	arcpy.AddField_management (outF6B, "LineID", "SHORT")
	arcpy.CalculateField_management(outF6B, "LineID", lineid,"PYTHON_9.3","")

	arcpy.AddField_management (outF6B, "WSTx", "DOUBLE")
	arcpy.CalculateField_management(outF6B, "WSTx", wstx,"PYTHON_9.3","")

	arcpy.AddField_management (outF6B, "WSTy", "DOUBLE")
	arcpy.CalculateField_management(outF6B, "WSTy", wsty,"PYTHON_9.3","")

	arcpy.AddField_management (outF6B, "WSCName", "TEXT")
	arcpy.CalculateField_management(outF6B, "WSCName", wscn,"PYTHON","")

	arcpy.AddField_management (outF6B, "ENTx", "DOUBLE")
	arcpy.CalculateField_management(outF6B, "ENTx", entx,"PYTHON","")

	arcpy.AddField_management (outF6B, "ENTy", "DOUBLE")
	arcpy.CalculateField_management(outF6B, "ENTy", enty,"PYTHON","")

	arcpy.AddField_management (outF6B, "ENCName", "TEXT")
	arcpy.CalculateField_management(outF6B, "ENCName", encn,"PYTHON","")

	# get vertices
	outF6D = outF6C.replace(".shp","2pt.shp")
	arcpy.FeatureVerticesToPoints_management (outF6C, outF6D, "ALL")
	arcpy.AddXY_management (outF6D)

	# using search cursor get all point coordinates
	SC = arcpy.SearchCursor(outF6D)
	xslvpts = []
	for row in SC:
		pt = [row.getValue("POINT_X"),row.getValue("POINT_Y")]
		xslvpts.append(pt)
	del row, SC,pt

	# using search cursor get all pseudo coordinates (plat, plon)
	SC = arcpy.SearchCursor(outF6B)
	pCords = []
	for row in SC:
		if setid == 16:
			pass							# # # UPDATE# # #
		elif setid == 17 or setid ==18:
			pass							# # # UPDATE# # #
		else:
			if float(row.getValue("pLon")) >= wstx and float(row.getValue("pLon")) <=entx:
				pt = [row.getValue("FID"),row.getValue("pLon"),row.getValue("pLat")]
				pCords.append(pt)
	del row, SC,pt

	# Create polyline shapefile(s) for linesegments  and their length
	out_path = os.path.dirname(outF6C)
	spRef = r"Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"
	# use variables setid,lineid,wstx,wsty,wscn,entx,enty,encn,pCords(FID,pLon,pLat),xslvpts(X,Y)
	if setid == 16:
		out_name1= DirCheck(out_path) + (os.path.basename(outF6C))[:3] + "_Dist2ECoast.shp"
		arcpy.CreateFeatureclass_management (out_path, out_name1, "POLYLINE", "", "DISABLED", "DISABLED", spRef)

								# # # UPDATE# # #

	elif setid == 17 or setid ==18:
		out_name1= DirCheck(out_path) + (os.path.basename(outF6C))[:3] + "_Dist2NCoast.shp"
		out_name2= DirCheck(out_path) + (os.path.basename(outF6C))[:3] + "_Dist2SCoast.shp"

								# # # UPDATE# # #

	else:
		out_name1= DirCheck(out_path) + (os.path.basename(outF6C))[:3] + "_Dist2ECoast.shp"
		out_name2= DirCheck(out_path) + (os.path.basename(outF6C))[:3] + "_Dist2WCoast.shp"
		entmp = DirCheck(out_path) + "entmp.csv"
		wstmp = DirCheck(out_path) + "wstmp.csv"
		entmpshplst = []
		wstmpshplst = []
		for pC in pCords:
			enf = []
			wsf = []
			for pt in xslvpts:
				if pC[1] < wstx and pC[1] < entx:
					if pt[0] < pC[1] or pt[0] > entx:
						pass
					elif pt >= pC[1] and pt[0] <= wstx:
						enf.append(pt)
						wsf.append(pt)
					elif pt[0] > pC[1] and pt[0] > wstx:
						enf.append(pt)
				elif pC[1] > wstx and pC[1] < entx:
					if pt[0] < wstx and pt[0] > entx:
						pass
					elif pt[0] <= pC[1]:
						wsf.append(pt)
					elif pt[0] >= pC[1]:
						enf.append(pt)
				elif pC[1] >= wstx and pC[1] >= entx:
					if pt[0] < wstx or pt[0] > pC[1]:
						pass
					elif pt[0] > wstx and pt[0] < entx:
						wsf.append(pt)
					elif pt[0] >= wstx and pt[0] >= entx:
						wsf.append(pt)
						enf.append(pt)
			spc = []
			spc.append(pC[1])
			spc.append(pC[2])
			enf.append(spc)
			wsf.append(spc)
			senf = []
			senf = sorted(enf, key=itemgetter(0))
			fw = open(entmp,'w')
			fw.write("lon,lat\n")
			for i in senf:
				wline = str(i[0]) + "," + str(i[1]) + "\n"
				fw.write(wline)
			fw.close()
			del fw,i,wline
			swsf = []
			swsf = sorted (wsf, key=itemgetter(0))
			fw = open(wstmp,'w')
			fw.write("lon,lat\n")
			for i in swsf:
				wline = str(i[0]) + "," + str(i[1]) + "\n"
				fw.write(wline)
			fw.close()
			del fw,i,wline

			# Make a feature layer
			try:
				t = t + 1
			except:
				t = 0
			enfl = "enFL" + str(t)
			wsfl = "wsFL" + str(t)
			arcpy.MakeXYEventLayer_management (entmp, "lon", "lat", enfl,spRef)
			arcpy.MakeXYEventLayer_management (wstmp, "lon", "lat", wsfl,spRef)
			# convert Points To Line
			outfolder = DirCheck(DirCheck(out_path) + "entmp")
			entmpofc = outfolder + str(pC[0]).zfill(3) +".shp"
			outfolder = DirCheck(DirCheck(out_path) + "wstmp")
			wstmpofc = outfolder + str(pC[0]).zfill(3) +".shp"
			entmpshplst.append(entmpofc)
			wstmpshplst.append(wstmpofc)
			arcpy.PointsToLine_management (enfl, entmpofc)
			arcpy.PointsToLine_management (wsfl, wstmpofc)
			arcpy.AddField_management (entmpofc, "ptFID", "SHORT")
			arcpy.AddField_management (wstmpofc, "ptFID", "SHORT")
			arcpy.CalculateField_management(entmpofc, "ptFID", pC[0],"PYTHON_9.3","")
			arcpy.CalculateField_management(wstmpofc, "ptFID", pC[0],"PYTHON_9.3","")
			arcpy.DeleteField_management (entmpofc, "ID")
			arcpy.DeleteField_management (wstmpofc, "ID")
			os.remove(entmp)
			os.remove(wstmp)
		print "No of files to merge for East Coas t= ", len(entmpshplst)
		arcpy.Merge_management(entmpshplst, out_name1)
		print "No of files to merge for West Coast = ", len(wstmpshplst)
		arcpy.Merge_management(wstmpshplst, out_name2)

# Project the shapefiles
pspRef = r"Coordinate Systems\Projected Coordinate Systems\World\WGS 1984 World Mercator.prj"
out_name11 = out_name1.replace(".shp", "_Proj.shp")
arcpy.Project_management (out_name1, out_name11, pspRef)
out_name21 = out_name1.replace(".shp", "_Proj.shp")
arcpy.Project_management (out_name2, out_name21, pspRef)

# Calculate length (Meters)
arcpy.AddField_management(out_name11, "Length","DOUBLE")
arcpy.CalculateField_management(out_name11, "Length", "!Shape!.length","PYTHON_9.3","")
arcpy.AddField_management(out_name21, "Length","DOUBLE")
arcpy.CalculateField_management(out_name21, "Length", "!Shape!.length","PYTHON_9.3","")
# join with Normals2XSectionLine_wET
# cleanup

"""-----------------------------------------------------------------------------
# End Time
-----------------------------------------------------------------------------"""
EndTime = time.asctime(time.localtime(time.time()))
print "\nEnd Time: ", EndTime,"\n"
