#-------------------------------------------------------------------------------
# Name:			SWMGenerator Part I.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 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)

""" Select input parameters """
inFC = BrowseFile("Select the Cross-Section Line Shapefile: ","C:/OES/","shp")
rl = "Total number of features in " + inFC + " = " + str(arcpy.GetCount_management (inFC)) + "\n"
fReport.write(rl)
print rl
del rl

""" Step 1: Duplicate Input File """
outF1 = DirCheck(inFC) + "Step1Out.shp"
arcpy.Copy_management(inFC, outF1)
rl = "Input file duplicated........"
fReport.write(rl)
print rl
del rl

""" Step 2: Add 'SetID' & 'LineID' fields and calculate LineID """
arcpy.AddField_management (outF1, "SetID", "SHORT")
# Get SetID from User & calculate field
##arcpy.CalculateField_management (outF1, "SetID", raw_input("Please enter the Box ID (interger)"))
arcpy.AddField_management (outF1, "LineID", "SHORT")
# Delelte the field "ID"
arcpy.DeleteField_management (outF1, "ID")
# Calculate distance of each feature to NprthPole
arcpy.Near_analysis (outF1, BrowseFile("Please browse for 'North Pole.shp'","C:/OES/","shp"))
# Delelte the field "NEAR_FID"
arcpy.DeleteField_management (outF1, "NEAR_FID")
# Make a list of Distance to NorthPole
DVals = ([row[0] for row in arcpy.da.SearchCursor(outF1, "NEAR_DIST")])
# Remove Duplicates if any
DVals = sets.Set(DVals)
# Sort and store in a new variable
sDVals = list(DVals)
sDVals.sort(reverse=False)
# Use UpdateCursor to calculate LineID
fields = ("LineID","NEAR_DIST")
with arcpy.da.UpdateCursor(outF1, fields) as cursor:
	for row in cursor:
##		print row
		row[0] = int(str(sDVals.index(row[1])+1).zfill(3))
		cursor.updateRow(row)
# Delelte the field "NEAR_DIST"
arcpy.DeleteField_management (outF1, "NEAR_DIST")
rl = "SetID and LineID fields added........"
fReport.write(rl)
print rl
del rl

""" Step 3: Add Terminal Points to the CS Lines """
# Intersect "World Major Continent Blocks.shp" and CSLines to get the intersecing points
in_features = [outF1,"C:/OES/Important Shapefiles/World Top50 Continental Blocks.shp"]
outF3A = DirCheck(inFC) + "Step3AOut.shp"	# ShapeFile generated (Point)
arcpy.Intersect_analysis (in_features, outF3A, "ALL", "", "POINT")

# Get FieldNames for this point shapefile and delete Extra FID fields
fns = []
fns = AllfFields(outF3A)
print fns, len(fns)
for items in fns:
	if "FID_" in items or "_FID" in items or "Area" in items:
		arcpy.DeleteField_management (outF3A, items)
	else:
		pass
del items

# Explode the multipart features
outF3B = DirCheck(inFC) + "Step3BOut.shp"	# ShapeFile generated (Point)
arcpy.MultipartToSinglepart_management (outF3A, outF3B)
arcpy.DeleteField_management (outF3B, "ORIG_FID")

# Add XY and rename fields
arcpy.AddXY_management (outF3B)
arcpy.AddField_management (outF3B, "lon", "DOUBLE")
expression = "!Point_X!"
arcpy.CalculateField_management (outF3B, "lon", expression, "PYTHON")
arcpy.DeleteField_management (outF3B, "Point_X")

arcpy.AddField_management (outF3B, "lat", "DOUBLE")
expression = "!Point_Y!"
arcpy.CalculateField_management (outF3B, "lat", expression, "PYTHON")
arcpy.DeleteField_management (outF3B, "Point_Y")

arcpy.AddField_management (outF3B, "CName", "TEXT")
expression = "!CONTINENT!"
arcpy.CalculateField_management (outF3B, "CName", expression, "PYTHON")
arcpy.DeleteField_management (outF3B, "CONTINENT")

rl = "lon, lat and CName (Continent) fields added........"
fReport.write(rl)
print rl
del rl

# Make a feature Layer
try:
	t = t + 1
except:
	t = 0
fl = "FL" + str(t)
arcpy.MakeFeatureLayer_management ( outF3B, fl)

# Make a list of SetID values
sids = [row[0] for row in arcpy.da.SearchCursor(outF3B, "SetID")]
sids = sets.Set(sids)
del row
for setid in sids:
##	print setid
	where_clause = "\"SetID\" = " + str(setid)
	##		print where_clause
	arcpy.SelectLayerByAttribute_management (fl, "NEW_SELECTION", where_clause)#Select features with certain SetID
	nosf = int(arcpy.GetCount_management(fl).getOutput(0))
	##		print "No of selected features: ", nosf
	linelist = [row[0] for row in arcpy.da.SearchCursor(fl, "LineID")]
	linelist = sets.Set(linelist)
	slinelist = []
	for items in linelist:
		slinelist.append(float(items))
	##		print slinelist #
	slinelist.sort()
	##		print slinelist #
	for slid in slinelist:
		where_clause = "\"SetID\" = " + str(setid) + " AND \"LineID\" = " + str(slid)#Select features with certain SetID and LineID
	##			print where_clause
		arcpy.SelectLayerByAttribute_management (fl, "NEW_SELECTION", where_clause)
		nosf = int(arcpy.GetCount_management(fl).getOutput(0))
	##			print "No of selected features: ", nosf
		cnamelist = [row[0] for row in arcpy.da.SearchCursor(fl, "CName")]
		cnamelist = sets.Set(cnamelist)
		for cname in cnamelist:
			where_clause = "\"SetID\" = " + str(setid) + " AND \"LineID\" = " + str(slid) + " AND \"CName\" = '" + cname + "'"#Select features with certain SetID, LineID and CName
	##				print where_clause
			arcpy.SelectLayerByAttribute_management (fl, "NEW_SELECTION", where_clause)
			nosf = int(arcpy.GetCount_management(fl).getOutput(0))
			if nosf >= 2:
	##					print "No of selected features: ", nosf
				# Decide to prepare for lonlist or latlist
				if setid == 17 or setid ==18:
##					print "Souther Ocean"
					latlist = [row[0] for row in arcpy.da.SearchCursor(fl, "lat")]
					latlist = sets.Set(latlist)
					slatlist = []
					for items in latlist:
						slatlist.append(float(items))
	##						print slatlist
					slatlist.sort()
	##						print slatlist
					if cname == "Australia":
						where_clause = "\"SetID\" = " + str(setid) + " AND \"LineID\" = " + str(slid) + " AND \"CName\" = '" + cname + "'" + " AND \"lat\" <> " + str(slatlist[0])
					elif cname == "Antarctica":
						where_clause = "\"SetID\" = " + str(setid) + " AND \"LineID\" = " + str(slid) + " AND \"CName\" = '" + cname + "'" + " AND \"lat\" <> " + str(slatlist[len(slatlist)-1])
##					print where_clause # # # #
				elif setid == 16:
##					print "Pacific"
					lonlist = [row[0] for row in arcpy.da.SearchCursor(fl, "lon")]
					lonlist = sets.Set(lonlist)
					slonlist = []
					for items in lonlist:
						slonlist.append(float(items))
	##						print slonlist
					slonlist.sort()
	##						print slonlist
					if cname == "South America":
						where_clause = "\"SetID\" = " + str(setid) + " AND \"LineID\" = " + str(slid) + " AND \"CName\" = '" + cname + "'" + " AND \"lon\" <> " + str(slonlist[0])
				else:
					lonlist = [row[0] for row in arcpy.da.SearchCursor(fl, "lon")]
					lonlist = sets.Set(lonlist)
					slonlist = []
					for items in lonlist:
						slonlist.append(float(items))
	##						print slonlist
					slonlist.sort()
	##						print slonlist
					if "America" in cname:
						where_clause = "\"SetID\" = " + str(setid) + " AND \"LineID\" = " + str(slid) + " AND \"CName\" = '" + cname + "'" + " AND \"lon\" <> " + str(slonlist[len(slonlist)-1])
					else:
						where_clause = "\"SetID\" = " + str(setid) + " AND \"LineID\" = " + str(slid) + " AND \"CName\" = '" + cname + "'" + " AND \"lon\" <> " + str(slonlist[0])
	##					print where_clause
				arcpy.SelectLayerByAttribute_management (fl, "NEW_SELECTION", where_clause)
				sftd = int(arcpy.GetCount_management(fl).getOutput(0))
##				print sftd, " features to delete"
				arcpy.DeleteFeatures_management (fl)
print outF3B, " cleaned"

""" Step 4: """
# Sort by "LineID"
outF4A = DirCheck(inFC) + "Step4AOut.shp"	# ShapeFile generated (Point)
arcpy.Sort_management(outF3B, outF4A, [["LineID", "ASCENDING"]])

# Export the attributes to a csv file
fieldNames = [X.name for X in arcpy.ListFields(outF4A)]

for fld in fieldNames:
	print "Field Name:  ",fld
	if "FID" in fld or fld =="Shape" or fld =="lon" or fld =="lat":
		fieldNames.remove(fld)
		print fld, " deleted..."

fieldNames = ["SetID","LineID","CName"]

outF4B = outF4A.replace("AOut.shp","BOut.csv")
# ExportXYv_stats (Input_Feature_Class, Value_Field, Delimiter, Output_ASCII_File, Add_Field_Names_to_Output)
arcpy.ExportXYv_stats(outF4A, fieldNames, "COMMA", outF4B, "ADD_FIELD_NAMES")

# Process the exported data for join with line features
f = open(outF4B, 'r')
fText = f.readlines()
f.close()

outF4C = outF4B.replace("BOut.csv","COut.csv")
fo = open(outF4C, 'w')
Header = "SetID,LineID,WSTx,WSTy,WSCName,ENTx,ENTy,ENCName\n"
##print Header

lNo = 1
WSTx = 0.0
ENTx = 0.0
WSTy = 0.0
ENTy = 0.0
WSCName = ""
ENCName = ""

for lines in fText:
	if lNo == 1:
		fo.write(Header)
	else:
		l = (lines.replace("\n","")).split(",")
		lon = l[0]
		lat = l[1]
		sid = l[2]
		lid = l[3]
		cname = l[4]
		if int(sid) == 16:
			wLine = sid + "," + lid + "," + str(-9999) + "," + str(-9999) + "," + "Not Applicable" + "," + lon + "," + lat + "," + cname + "\n"
##			print wline
			fo.write(wLine)
		elif int(sid) == 17 or int(sid) == 18:
			if cname == "Australia":
				sid1 = sid
				lid1 = lid
				ENTx = lon
				ENTy = lat
				ENCName = cname
			elif cname == "Antarctica":
				sid2 = sid
				lid2 = lid
				WSTx = lon
				WSTy = lat
				WSCName = cname
		else:
			if "America" in cname:
				sid1 = sid
				lid1 = lid
				WSTx = lon
				WSTy = lat
				WSCName = cname
			elif "Europe" in cname or "Africa" in cname:
				sid2 = sid
				lid2 = lid
				ENTx = lon
				ENTy = lat
				ENCName = cname
		if ENTx != 0.0 and ENTy != 0.0 and ENCName != "" and WSTx != 0.0 and WSTy != 0.0 and WSCName != "":
			if sid1 == sid2 and lid1 == lid2:
				wline = sid + "," + lid + "," + WSTx + "," + WSTy + "," + WSCName + "," + ENTx + "," + ENTy + "," + ENCName + "\n"
##				print wline
				fo.write(wline)
				WSTx = 0.0
				ENTx = 0.0
				WSTy = 0.0
				ENTy = 0.0
				WSCName = ""
				ENCName = ""
	lNo = lNo + 1
fo.close()

""" Step 5: Joining the endpoints to cross section lines"""
# Duplicate Step1Out and join to Step4C table using LineID
outF5A = DirCheck(inFC) + "Step5AOut.shp"
arcpy.Copy_management(outF1, outF5A)
rl = outF1 + " duplicated........"
fReport.write(rl)
print rl
del rl

# Join this outF5A with out4c csv
try:
	t = t + 1
except:
	t = 0
fl = "FL" + str(t)
arcpy.MakeFeatureLayer_management (outF5A, fl)
arcpy.AddJoin_management (fl, "LineID", outF4C, "LineID", "KEEP_ALL")

# export joined shapefile
outF5B = DirCheck(inFC) + "Step5BOut.shp"	# ShapeFile generated (Polyline)
arcpy.CopyFeatures_management( fl, outF5B)

# remove join
arcpy.RemoveJoin_management (fl)
del fl

# add fileds
arcpy.AddField_management (outF5B, "SetID", "SHORT")
arcpy.AddField_management (outF5B, "LineID", "SHORT")
arcpy.AddField_management (outF5B, "WSTx", "DOUBLE")
arcpy.AddField_management (outF5B, "WSTy", "DOUBLE")
arcpy.AddField_management (outF5B, "WSCName", "TEXT")
arcpy.AddField_management (outF5B, "ENTx", "DOUBLE")
arcpy.AddField_management (outF5B, "ENTy", "DOUBLE")
arcpy.AddField_management (outF5B, "ENCName", "TEXT")
# Calculate field values
expression = "!Step4COut_!"
arcpy.CalculateField_management (outF5B, "SetID", expression, "PYTHON")
expression = "!Step4COut1!"
arcpy.CalculateField_management (outF5B, "LineID", expression, "PYTHON")
expression = "!Step4COu_1!"
arcpy.CalculateField_management (outF5B, "WSTx", expression, "PYTHON")
expression = "!Step4COu_2!"
arcpy.CalculateField_management (outF5B, "WSTy", expression, "PYTHON")
expression = "!Step4COu_3!"
arcpy.CalculateField_management (outF5B, "WSCName", expression, "PYTHON")
expression = "!Step4COu_4!"
arcpy.CalculateField_management (outF5B, "ENTx", expression, "PYTHON")
expression = "!Step4COu_5!"
arcpy.CalculateField_management (outF5B, "ENTy", expression, "PYTHON")
expression = "!Step4COu_6!"
arcpy.CalculateField_management (outF5B, "ENCName", expression, "PYTHON")
# Delete fields
arcpy.DeleteField_management (outF5B, "Step5AOut_")
arcpy.DeleteField_management (outF5B, "Step5AOut1")
arcpy.DeleteField_management (outF5B, "Step4COut_")
arcpy.DeleteField_management (outF5B, "Step4COut1")
arcpy.DeleteField_management (outF5B, "Step4COu_1")
arcpy.DeleteField_management (outF5B, "Step4COu_2")
arcpy.DeleteField_management (outF5B, "Step4COu_3")
arcpy.DeleteField_management (outF5B, "Step4COu_4")
arcpy.DeleteField_management (outF5B, "Step4COu_5")
arcpy.DeleteField_management (outF5B, "Step4COu_6")

""" Section 6: dicing up for cross section """
outF5B = "C:/OES/Important Shapefiles/Step5BOut.shp"
# Select CSline by LineID
try:
	t = t + 1
except:
	t = 0
fl1 = "FL" + str(t)
arcpy.MakeFeatureLayer_management (outF5B, fl1)
try:
	t = t + 1
except:
	t = 0
fl2 = "FL" + str(t)
arcpy.MakeFeatureLayer_management ("C:/OES/PaleoBathymetry/eTopo/eTopo4XSentions.shp", fl2)
i = 1
while i<316:
	outF6A = DirCheck(inFC) + "diced/" + str(i) + "_XSectionLine.shp"
	where_clause = "\"LineID\" = " + str(i)
	arcpy.SelectLayerByAttribute_management (fl1, "NEW_SELECTION", where_clause)
	arcpy.CopyFeatures_management(fl1, outF6A)
	outF6B = DirCheck(inFC) + "diced/" + str(i) + "_LineETPoints.shp"
	arcpy.SelectLayerByLocation_management (fl2, "WITHIN_A_DISTANCE", fl1, "0.25", "NEW_SELECTION")
	arcpy.CopyFeatures_management(fl2, outF6B)
	out_path = os.path.dirname(outF6B)
	out_name = str(i) + "_Normals2XSectionLine.shp"
	spRef = r"Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"
	arcpy.CreateFeatureclass_management (out_path, out_name, "POLYLINE", "", "DISABLED", "DISABLED", spRef)
	print "LineID ", str(i), " processed......"
	i = i + 1

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