#####
##	BEGIN:	Set general configuration options
##	Configuration:
root_folder = r"C:\Users\preed\Desktop\Thesis\Comparisons\Scripts\Create Courser Census Data"

base_census_name = "Haiti_2015_population_estimate_adm3"

##	This number is set so only a fraction of the original total is
##		removed before running dissolve.  A higher number will decrease
##		speed, but with the advantage of not leaving behind isolated units that
##		can't be merged because all of their neighbors are already involved
##		in a dissolution:
reduction_factor = 10

##	Random number seed:
seed = 2002

##	Target admin unit count:
target_count = 380

##	Sort and merge, smallest first if True, or by random if False:
merge_smallest = True

##	Set to true to create a fully-dissolved output layer:
merge_all = True

##	END:	Set general configuration options
#####



#####

##	BEGIN: Import packages and set GeoProcessing environment

import arcpy
import os

##	Should we overwrite any already existing derived data:
overwrite = True
arcpy.env.overwriteOutput = overwrite

##	END: Import packages and set GeoProcessing environment
#####



#####
##	BEGIN: Define utility functions

##	Create a directory function to check for and create data directories if they don't exist:
def ensure_dir(d):
	if not os.path.isdir(d):
		os.makedirs(d)
	return d

##	END: Define utility functions
#####



#####
##	BEGIN: Data processing

##	Create base folder and data names from above information:
input_data = root_folder + "/" + base_census_name + ".shp"
output_data_base = root_folder + "/" + base_census_name

##	Copy our original dataset to our output folder:
eliminate_data = output_data_base + "_eliminate.shp"
arcpy.Copy_management(input_data, eliminate_data, "")
input_data = eliminate_data


##	Remove un-needed fields:
fields = arcpy.ListFields(input_data)

required_fields = ["FID", "Shape", "ADMINID", "YEARPOP", "ADMINPOP", "Area", "Random", "Select", "Dissolve"]
delete_fields = []
for field in fields:
	if field.name not in required_fields:
		delete_fields.append(field.name)

arcpy.DeleteField_management(input_data, delete_fields)


##	Add a fields for random sorting and selection to our shapefile:
arcpy.AddField_management(input_data,"Area","FLOAT","#","#","#","#","NON_NULLABLE","NON_REQUIRED","#")
arcpy.AddField_management(input_data,"Random","FLOAT","#","#","#","#","NON_NULLABLE","NON_REQUIRED","#")
arcpy.AddField_management(input_data,"Select","SHORT","#","#","#","#","NON_NULLABLE","NON_REQUIRED","#")
arcpy.AddField_management(input_data,"Dissolve","LONG","#","#","#","#","NON_NULLABLE","NON_REQUIRED","#")


##	We can use the shape area to sort our aggregation, dissolving
##		smallest units first (or not depending on how the sort is set in our
##		cursors below:
arcpy.CalculateField_management(input_data,"Area","!SHAPE.AREA!","PYTHON","")


original_count = int(arcpy.GetCount_management(input_data)[0])
new_count = original_count
old_count = 0
i = 0
j = 0


while new_count != target_count and i != -1:
	print("Running elimination on " + str(new_count) + " census units...")

	##	Set default dissolve field to be equal to the remaining FIDs:
	arcpy.CalculateField_management(input_data,"Dissolve",'!FID!',"PYTHON","")

	##	Set the Select field to -1.  -1 will indicate non-selected, viable
	##		focal polgyon for dissolution, -2 will indicate that this polygon
	##		is absorbing another polygon and should not itself removed, and
	##		any value 0 or greater will indicate the polgyon that is being
	##		dissolved into:
	arcpy.CalculateField_management(input_data,"Select",'-1',"PYTHON","")

	##	Generate a set of random numbers to sort by if we are
	##		eliminating rows at random:
	arcpy.CalculateField_management(input_data,"Random","set_random()","PYTHON","import random\nrandom.seed(" + str(seed + new_count) + ")\ndef set_random():\n\treturn random.random()")


	##	ArcGIS 10.0 does not accept keyword arguments, so you have to specify
	##		them all...
	#UpdateCursor (dataset, {where_clause}, {spatial_reference}, {fields}, {sort_fields})
	if merge_smallest:
		focal_rows = arcpy.SearchCursor(input_data, '"Select" = -1', "", "","Area A")
	else:
		focal_rows = arcpy.SearchCursor(input_data, '"Select" = -1', "", "","Random A")


	# Make a layer from the feature class
	arcpy.MakeFeatureLayer_management(input_data, "lyr")

	i = 0
	j = 0
	while i in range(0, (original_count - target_count) / reduction_factor) and (new_count - i > target_count):
		print("Processing polygon: " + str(j))
		focal_row = focal_rows.next()

		if (focal_row != None):
			arcpy.SelectLayerByAttribute_management("lyr","NEW_SELECTION",'"FID"=' + str(focal_row.getValue("FID")))

			##	Select the surrounding units from our current, focal sampled unit:
			arcpy.SelectLayerByLocation_management("lyr","BOUNDARY_TOUCHES","lyr","#","NEW_SELECTION")

			##	Check to see if we have any touching features:
			focal_fid = focal_row.getValue("FID")
			merged_fid = None

			if (int(arcpy.GetCount_management("lyr").getOutput(0)) == 1):
				##	If we are an island, then calculate nearest polygon:
				print("Calculating nearest census units...")
				arcpy.Near_analysis("lyr", input_data)

				island_rows = arcpy.SearchCursor("lyr", '"FID" = ' + str(focal_fid), "", "","")
				island_row = island_rows.next()

				##	First we check to see if we have touching polygons and if not
				##		choose the nearest polygon to dissolve with:
				if (island_row.getValue("NEAR_DIST") > 0):
					merged_fid = island_row.getValue("NEAR_FID")
				del island_row
				del island_rows

			else:
				##	Our nearest polygon shares a border so let's determine which ones
				##		touch and choose the longest shared border to merge with:
				##	Select the focal unit:

				##	Remove any from the selection that are currently already being
				##		dissolved:
				arcpy.SelectLayerByAttribute_management("lyr","REMOVE_FROM_SELECTION",'"FID"=' + str(focal_row.getValue("FID")) + ' OR "Select" <> -1')

				##	Update the focal polygon's Dissolve field for the first touching
				##		polygon, sorted by smallest to largest:
				merge_rows = arcpy.SearchCursor("lyr", "", "", "", "Area A")
				merge_row = merge_rows.next()


				if (merge_row != None):

					##	Search through touching polygons and find the polygon with the
					##		longest shared border:

					##	Get the geometry of our focal unit:
					focal_shape = focal_row.getValue("SHAPE")

					focal_fid = focal_row.getValue("FID")
					merged_fid = None
					longest_border_len = 0

					while merge_row != None:
						merge_shape = merge_row.getValue("SHAPE")
						##	Calculate the intersection of the two geometries along two
						##		dimensions (the 2) indicating we want a linear output:
						##	NOTE: This will not work pre-10.1!!!
						intersection = focal_shape.intersect(merge_shape, 2)
						if intersection.getLength() > longest_border_len:
							merged_fid = merge_row.getValue("FID")
							longest_border_len = intersection.getLength()
						merge_row = merge_rows.next()

				del merge_row
				del merge_rows

			if merged_fid != None:
				##	Set the focal polygon's Select and Dissolve field:
				update_rows = arcpy.UpdateCursor(input_data, '"FID"='+str(focal_row.getValue("FID")), "", "","")
				update_row = update_rows.next()
				update_row.setValue("Select", merged_fid)
				update_row.setValue("Dissolve", merged_fid)
				update_rows.updateRow(update_row)
				del update_row
				del update_rows

				##	Set the polygon to be dissolved into's Select field so it's not
				##		dissolved out later on:
				update_rows = arcpy.UpdateCursor(input_data, '"FID"='+str(merged_fid), "", "","")
				update_row = update_rows.next()
				update_row.setValue("Select", -2)
				update_rows.updateRow(update_row)
				del update_row
				del update_rows

				print("Reassigned polygon " + str(i) + ":  " + str(focal_row.getValue("FID")) + " to " + str(merged_fid))
				i += 1

			else:
				print("WARNING: Polygon has no remaining neighbors to merge with but is not an island...  Skipped!")

								##	Set the focal polygon's Select:
				update_rows = arcpy.UpdateCursor(input_data, '"FID"='+str(focal_row.getValue("FID")), "", "","")
				update_row = update_rows.next()
				update_row.setValue("Select", -3)
				update_rows.updateRow(update_row)
				del update_row
				del update_rows

			##	Last, recreate our SearchCursor so the updated data is included:
			del focal_row
			del focal_rows

			if merge_smallest:
				focal_rows = arcpy.SearchCursor(input_data, '"Select" = -1', "", "","Area A")
			else:
				focal_rows = arcpy.SearchCursor(input_data, '"Select" = -1', "", "","Random A")

		else:
			##	If we are out of rows then stop:
			if focal_row == None or i == original_count - target_count:
				new_count = new_count - i - 1
				i = -1

		j += 1

	del focal_rows

	output_data = output_data_base + "_" + str(new_count-i) + ".shp"
	arcpy.Dissolve_management(input_data, output_data,"Dissolve","ADMINID FIRST;YEARPOP FIRST;ADMINPOP SUM;Area SUM;Random FIRST;Select MAX","MULTI_PART","DISSOLVE_LINES")

	##	We have the irritating problem that the Dissolve tool renames fields
	##		according to their merge function *and* ArcGIS has no way of
	##		renaming fields, so have to create a new field and copy values:
	fields = arcpy.ListFields(output_data)

	old_fields = ["FIRST_ADMI", "FIRST_YEAR", "SUM_ADMINP", "SUM_Area", "FIRST_Rand", "MAX_Select", "FIRST_Diss"]
	new_fields = ["ADMINID", "YEARPOP", "ADMINPOP", "Area", "Random", "Select"]
	new_fields_type = ["LONG", "LONG", "LONG", "FLOAT", "FLOAT", "SHORT"]

	for i in range(0,len(new_fields)):
		arcpy.AddField_management(output_data,new_fields[i],new_fields_type[i],"#","#","#","#","NON_NULLABLE","NON_REQUIRED","#")
		arcpy.CalculateField_management(output_data,new_fields[i],"!" + old_fields[i] + "!","PYTHON","")

	arcpy.DeleteField_management(output_data, old_fields)


	##	Calculate our new count:
	input_data = output_data
	new_count = int(arcpy.GetCount_management(input_data)[0])


##	Finally, do one last dissolve to create a single polygon result:
if merge_all:
	output_data = output_data_base + "_1.shp"
	arcpy.Dissolve_management(input_data, output_data,"","ADMINID FIRST;YEARPOP FIRST;ADMINPOP SUM;Area SUM;Random FIRST;Select MAX","MULTI_PART","DISSOLVE_LINES")

	##	We have the irritating problem that the Dissolve tool renames fields
	##		according to their merge function *and* ArcGIS has no way of
	##		renaming fields, so have to create a new field and copy values:
	fields = arcpy.ListFields(output_data)

	old_fields = ["FIRST_ADMI", "FIRST_YEAR", "SUM_ADMINP", "SUM_Area", "FIRST_Rand", "MAX_Select", "FIRST_Diss"]
	new_fields = ["ADMINID", "YEARPOP", "ADMINPOP", "Area", "Random", "Select"]
	new_fields_type = ["LONG", "LONG", "LONG", "FLOAT", "FLOAT", "SHORT"]

	for i in range(0,len(new_fields)):
		arcpy.AddField_management(output_data,new_fields[i],new_fields_type[i],"#","#","#","#","NON_NULLABLE","NON_REQUIRED","#")
		arcpy.CalculateField_management(output_data,new_fields[i],"!" + old_fields[i] + "!","PYTHON","")

	arcpy.DeleteField_management(output_data, old_fields)


##	Remove our eliminate temporary data shapefile:
arcpy.Delete_management(eliminate_data)


##	END: Data processing
#####

