polygon_divider.py

# -*- coding: utf-8 -*-
"""
/***************************************************************************
 PolygonDivider
								 A QGIS plugin
 Divides polygons into smaller 'squareish' polygons of a specified size
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
							  -------------------
		begin                : 2020-02-24
		git sha              : $Format:%H$
		copyright            : (C) 2023 Lune Geographic
		email                : info@lunegeographic.co.uk
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

 * This script divides a polygon into squareish sections of a specified size
 *
 * ERROR WE CAN FIX BY ADJUSTING TOLERANCE / N_SUBDIVISIONS:
 *  - Bracket is smaller than tolerance: the shape got smaller? OR got dramatically bigger. Can we check this? This is where we just want to cut at the last location that worked and re-calculate the division stuff.
 *
 * TODO'S:
 *  - Needs some form of form validation
 *  - There are two deprecated functions in this code
 *  - Pull out repeated file writing stuff into a function
 *  - Where / how often should we calculate the desired area?
 *  - How should we be dealing with reversing direction for subdivision?
 *  - Need to re-think how to deal with problems in subdivision - perhaps need to calculate all subdivisions then save all at once so we can roll back?
 *  - Look at saving last good bounds to narrow search interval after adjusting tolerance? Maybe use bisection to minimise the adjustment in tolerance?
 *  - This could be run in multiple parallel tasks (all strips calculated then divided up into processes to get squareish sections from them)
 *  - More minor TODO's throughout the code...
 *
 * @author jonnyhuck
 *
"""

import sys
import os.path
from uuid import uuid4
from .resources import *
from qgis.utils import iface
from qgis.PyQt.QtGui import QIcon
from qgis.PyQt.QtWidgets import QAction, QFileDialog
from .polygon_divider_dialog import PolygonDividerDialog
from qgis.PyQt.QtCore import QVariant, QSettings, QTranslator, QCoreApplication
from qgis.core import Qgis, QgsGeometry, QgsPointXY, QgsPoint, QgsField, QgsTask, QgsFeature, QgsVectorLayer, \
	QgsVectorFileWriter, QgsProject, QgsMessageLog, QgsApplication, QgsWkbTypes


class BrentError(Exception):
	"""
	* Simple class used for exceptions from Brent's Method.
	"""
	def __init__(self, value):
		self.value = value
	def __str__(self):
		return repr(self.value)


# Global label for task messages
MESSAGE_CATEGORY = 'POLYGON DIVIDER'


class PolygonDividerTask(QgsTask):
	"""
	* Class extending QgsTask to handle on different thread
	"""

	''' INHERITED CLASS FUNCTIONS '''

	def __init__(self, layer, outFilePath, target_area, absorb_flag, direction, tolerance):
		"""
		* Initialise the thread
		"""
		# call constructor for QsgThread
		super().__init__("Polygon Division", QgsTask.CanCancel)

		# load arguments into class variables
		self.layer = layer
		self.outFilePath = outFilePath
		self.target_area = target_area
		self.absorb_flag = absorb_flag
		self.direction = direction
		self.tolerance = tolerance
		self.exception = None


	def finished(self, result):
		"""
		* This function is automatically called when the task has completed (successfully or not).
		* It called from the main thread, so it's safe to do GUI operations and raise Python exceptions here.
		* result is the return value from self.run.
		"""
		# success
		if result:

			# Notify User
			QgsMessageLog.logMessage('Polygon Division completed', MESSAGE_CATEGORY, Qgis.Success)
			iface.messageBar().pushMessage("Success!", 'Polygon Division completed', level=Qgis.Success, duration=3)

			# finally, open the resulting file and return it
			layer = iface.addVectorLayer(self.outFilePath, 'Divided Polygon', "ogr")

			# alert if invalid
			if not layer.isValid():
				raise Exception("Output Dataset Invalid")

		# failure
		else:
			# failed without exception (user cancel)
			if self.exception is None:
				QgsMessageLog.logMessage('Polygon Division exited without exception', MESSAGE_CATEGORY, Qgis.Warning)
				iface.messageBar().pushMessage("Warning:", 'Polygon Division Cancelled', level=Qgis.Warning, duration=3)

			# failed with exception (raise exception)
			else:
				QgsMessageLog.logMessage(f'Polygon Division Exception: {self.exception}', MESSAGE_CATEGORY, Qgis.Critical)
				iface.messageBar().pushMessage("Error:", f'Polygon Division Failed: {str(self.exception)}', level=Qgis.Critical)
				# raise self.exception


	def cancel(self):
		"""
		* This function is called when the task is cancelled
		"""

		# log that the task was cancelled
		QgsMessageLog.logMessage("Polygon Division was manually cancelled by the user", MESSAGE_CATEGORY, Qgis.Info)

		# run the cancel method from QgsTask
		super().cancel()


	''' OPERATIONAL FUNCTIONS '''


	def brent(self, xa, xb, xtol, ftol, max_iter, geom, fixedCoord1, fixedCoord2, targetArea, horizontal, forward):
		"""
		* Brent's Method, following Wikipedia's article algorithm.
		*
		* - xa is the lower bracket of the interval of the solution we search.
		* - xb is the upper bracket of the interval of the solution we search.
		* - xtol is the minimum	 permitted width (in map units) of the interval before we give up.
		* - ftol is the required precision of the solution.
		* - max_iter is the maximum allowed number of iterations.
		* - geom is the geometry we are dividing.
		* - fixedCoord1 is the the first coordinate in the fixed dimension.
		* - fixedCoord2 is the the second coordinate in the fixed dimension.
		* - targetArea is the desired area of the section to cut off geom.
		* - horizontal or vertical cut - True / False respectively
		* - forward (left-right or bottom top) cut or backward (right-left or top-bottom) - True / False respectively
		"""

		''' SETUP '''

		# standard for iterative algorithms
		EPS = sys.float_info.epsilon

		''' BASIC ERROR CHECKING (INTERVAL VALIDITY) '''

		# check that the bracket's interval is sufficiently big for this computer to work with.
		if abs(xb - xa) < EPS:
			# raise BrentError("Initial bracket smaller than system epsilon.")
			self.exception = BrentError("Initial bracket smaller than system epsilon.")
			return False

		# check lower bound
		fa = self.f(xa, geom, fixedCoord1, fixedCoord2, targetArea, horizontal, forward)   # first function call
		if abs(fa) < ftol:
			# raise BrentError("Root is equal to the lower bracket")
			self.exception = BrentError("Root is equal to the lower bracket")
			return False

		# check upper bound
		fb = self.f(xb, geom, fixedCoord1, fixedCoord2, targetArea, horizontal, forward)   # second function call
		if abs(fb) < ftol:
			# raise BrentError("Root is equal to the upper bracket")
			self.exception = BrentError("Root is equal to the upper bracket")
			return False

		# check if the root is bracketed.
		if fa * fb > 0.0:	# this is checking for different signs (to be sure we are either side of 0)
			# raise BrentError("Root is not bracketed.")
			self.exception = BrentError("Root is not bracketed.")
			return False

		''' START CALCULATION '''

		# if the area from a is smaller than b, switch the values
		if abs(fa) < abs(fb):
			xa, xb = xb, xa
			fa, fb = fb, fa

		# initialise c at a (therefore at the one with the biggest area to the right of it)
		xc, fc = xa, fa

		# init mflag
		mflag = True

		# do until max iterations is reached
		for i in range(max_iter):

			# try to calculate `xs` by using inverse quadratic interpolation...
			if fa != fc and fb != fc:
				xs = (xa * fb * fc / ((fa - fb) * (fa - fc)) + xb * fa * fc / ((fb - fa) * (fb - fc)) + xc * fa * fb / ((fc - fa) * (fc - fb)))
			else:
				# ...if you can't, use the secant rule.
				xs = xb - fb * (xb - xa) / (fb - fa)

			# check if the value of `xs` is acceptable, if it isn't use bisection.
			if ((xs < ((3 * xa + xb) / 4) or xs > xb) or
				(mflag == True	and (abs(xs - xb)) >= (abs(xb - xc) / 2)) or
				(mflag == False and (abs(xs - xb)) >= (abs(xc - d) / 2)) or
				(mflag == True	and (abs(xb - xc)) < EPS) or
				(mflag == False and (abs(xc - d)) < EPS)):

				# overwrite unacceptable xs value with result from bisection
				xs = (xa + xb) / 2
				mflag = True
			else:
				mflag = False

			''' THE ABOVE BLOCK USED BRENT'S METHOD TO GET A SUGGESTED VALUE FOR S,
				THE BELOW BLOCK CHECKS IF IT IS GOOD, AND IF NOT SEEKS A NEW VALUE '''

			# get the value from f using the new xs value
			fs = self.f(xs, geom, fixedCoord1, fixedCoord2, targetArea, horizontal, forward) # repeated function call

			# if the value (ideally 0) is less than the specified tolerance, return
			if abs(fs) < ftol:
				return xs

			# if the bracket has become smaller than the tolerance (but the value wasn't reached, something is wrong)
			# this can indicate the 'W' condition, where decreasing the interval increases the size of the resulting area
			if abs(xb - xa) < xtol:
				# raise BrentError("Bracket is smaller than tolerance.")
				self.exception = BrentError("Bracket is smaller than tolerance.")
				return False

			# d is assigned for the first time here; it won't be used above on the first iteration because mflag is set
			d = xc	# it is just used in Brent's checks, not in the calculation per se

			# move c to b
			xc, fc = xb, fb

			# move one of the interval edges to the new point, such that zero remains within the interval
			# if the areas from a and s (current result) are same sign, move b to s, otherwise, move a to s
			if fa * fs < 0:			# different signs
				xb, fb = xs, fs
			else:					# same sign
				xa, fa = xs, fs

			# if the area from a is smaller than b, switch the values
			if abs(fa) < abs(fb):
				xa, xb = xb, xa
				fa, fb = fb, fa

		# this isn't as good (ran out of iterations), but seems generally fine
		return xs	# NB: increasing the number of iterations doesn't seem to get any closer


	def splitPoly(self, polygon, splitter, horizontal, forward):
		"""
		* Split a Polygon with a LineString
		* Returns exactly two polygons notionally referred to as being 'left' and 'right' of the cutline.
		* The relationship between them is that the 'right' polygon (the chip) is notionally cut from the 'left' one (the potato).
		"""

		# need to take a deep copy for the incoming polygon, as splitGeometry edits it directly...
		poly = QgsGeometry(polygon)

		'''
		TODO: SPLITGEOMETRY IS DEPRECATED, but...

		TypeError: index 0 has type 'QgsPoint' but 'QgsPointXY' is expected 
		Traceback (most recent call last):
		File "/Users/jonnyhuck/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/Submission/polygon_divider.py", line 126, in finished
			raise self.exception
		File "/Users/jonnyhuck/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/Submission/polygon_divider.py", line 611, in run
			sqArea = self.getSliceArea(interval[1] - sq + buffer, poly, fixedCoords[0], fixedCoords[1], horizontal_flag, forward_flag)	# cutting from top/right
		File "/Users/jonnyhuck/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/Submission/polygon_divider.py", line 429, in getSliceArea
			left, right, residual = self.splitPoly(poly, splitter, horizontal, forward)
		File "/Users/jonnyhuck/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/Submission/polygon_divider.py", line 282, in splitPoly
			res, polys, topolist = poly.splitGeometry(splitter, False)
		TypeError: index 0 has type 'QgsPoint' but 'QgsPointXY' is expected

		This is due to...
			https://github.com/qgis/QGIS/issues/37821

		See deprecation notes here:
			https://api.qgis.org/api/deprecated.html
		'''
		# parse splitter into required objects
		# splitter = [QgsPointXY(splitter[0][0], splitter[0][1]), QgsPointXY(splitter[1][0], splitter[1][1])]
		splitter = [QgsPoint(splitter[0][0], splitter[0][1]), QgsPoint(splitter[1][0], splitter[1][1])]

		# split poly (polygon) by splitter (line)
		# http://gis.stackexchange.com/questions/114414/cannot-split-a-line-using-qgsgeometry-splitgeometry-in-qgis
		res, polys, topolist = poly.splitGeometry(splitter, False)

		# add poly (which might be a multipolygon) to the polys array
		if poly.isMultipart():
			## TODO: I think that this is where we might be getting an odd error on no absorb
			# if the feature is a multipolygon, explode into separate polygons for individual processing
			multiGeom = poly.asMultiPolygon()
			for i in multiGeom:
				polys.append(QgsGeometry().fromPolygonXY(i))
		else:
			# ...OR load the feature into a list of one (it may be extended in the course of splitting if we create noncontiguous offcuts) and loop through it
			polys.append(poly)

		''' SORT RIGHT, LEFT, RESIDUAL '''

		# verify that it worked and that more than one polygon was returned
		if res == 0 and len(polys) >1:
			if forward: ### from bottom left
				if horizontal:	## cut from the bottom

					# left is the top one
					maxy = float('-inf')
					for i in range(len(polys)):
						p = polys[i].boundingBox().yMaximum()
						if p > maxy:
							maxy = p
							maxyi = i
					left = polys.pop(maxyi)

					# right is the bottom one
					miny = float('inf')
					for i in range(len(polys)):
						p = polys[i].boundingBox().yMinimum()
						if p < miny:
							miny = p
							minyi = i
						elif p == miny:		# if there is a tie for which is the rightest, get the rightest in the other dimension
							if polys[i].boundingBox().xMinimum() < polys[minyi].boundingBox().xMinimum():	# left
								minyi = i
					right = polys.pop(minyi)

				else:	## cutting from the left

					# left is the rightest one
					maxx = float('-inf')
					for i in range(len(polys)):
						p = polys[i].boundingBox().xMaximum()
						if p > maxx:
							maxx = p
							maxxi = i
					left = polys.pop(maxxi)

					# right is the leftest one
					minx = float('inf')
					for i in range(len(polys)):
						p = polys[i].boundingBox().xMinimum()
						if p < minx:
							minx = p
							minxi = i
						elif p == minx:		# if there is a tie for which is the rightest, get the rightest in the other dimension
							if polys[i].boundingBox().yMinimum() < polys[minxi].boundingBox().yMinimum():	# bottom
								minxi = i
					right = polys.pop(minxi)

			else:	### cut from top / right (forward_flag is false)

				if horizontal:	## cut from the top

					# left is the bottom one
					miny = float('inf')
					for i in range(len(polys)):
						p = polys[i].boundingBox().yMinimum()
						if p < miny:
							miny = p
							minyi = i
					left = polys.pop(minyi)

					# right is the top one
					maxy = float('-inf')
					for i in range(len(polys)):
						p = polys[i].boundingBox().yMaximum()
						if p > maxy:
							maxy = p
							maxyi = i
						elif p == maxy:		# if there is a tie for which is the rightest, get the rightest in the other dimension
							if polys[i].boundingBox().xMaximum() > polys[maxyi].boundingBox().xMaximum():
								maxyi = i
					right = polys.pop(maxyi)

				else:	## cutting from the right

					# left is the leftest one
					minx = float('inf')
					for i in range(len(polys)):
						p = polys[i].boundingBox().xMinimum()
						if p < minx:
							minx = p
							minxi = i
					left = polys.pop(minxi)

					# right is the rightest one
					maxx = float('-inf')
					for i in range(len(polys)):
						p = polys[i].boundingBox().xMaximum()
						if p > maxx:
							maxx = p
							maxxi = i
						elif p == maxx:		# if there is a tie for which is the rightest, get the rightest in the other dimension
							if polys[i].boundingBox().yMaximum() > polys[maxxi].boundingBox().yMaximum():
								maxxi = i
					right = polys.pop(maxxi)

			# work out if any remaining polygons are contiguous with left or not
			contiguous = []
			noncontiguous = []
			if len(polys) > 0:
				for j in polys:
					if left.touches(j):
						contiguous.append(j)
					else:
						noncontiguous.append(j)

				# join all contiguous parts back to left
				if len(contiguous) > 0:
					contiguous += [left]
					left = QgsGeometry.unaryUnion(contiguous)

			# return the two sections (left is the potato, right is the chip...), plus any noncontiguous polygons
			return left, right, noncontiguous
		else:
			# log error
			QgsMessageLog.logMessage("FAIL: Polygon division failed.", level=Qgis.Critical)


	def getSliceArea(self,sliceCoord, poly, fixedCoord1, fixedCoord2, horizontal, forward):
		"""
		* Splits a polygon to a defined distance from the minimum value for the selected dimension and
		*  returns the area of the resultng polygon
		"""

		# construct a list of points representing a line by which to split the polygon
		if horizontal:
			splitter = [(fixedCoord1, sliceCoord), (fixedCoord2, sliceCoord)] # horizontal split
		else:
			splitter = [(sliceCoord, fixedCoord1), (sliceCoord, fixedCoord2)] # vertical split

		# split the polygon
		left, right, residual = self.splitPoly(poly, splitter, horizontal, forward)

		# return the area of the bit you cut off
		return right.area()


	def f(self,sliceCoord, poly, fixedCoord1, fixedCoord2, targetArea, horizontal, forward):
		"""
		* Split a Polygon with a LineString, returning the area of the polygon to the right of the split
		* returns the area of the polygon on the right of the splitter line
		"""

		# return the difference between the resulting polygon area (right of the line) and the desired area
		return self.getSliceArea(sliceCoord, poly, fixedCoord1, fixedCoord2, horizontal, forward) - targetArea


	''' THE MAIN FUNCTION FOR THE TASK'''


	def run(self):
		"""
		* Actually do the processing
		"""

		# This whole function is in a catch-all try statement to avoid QGIS burying exceptions.
		try:

			# setup for progress bar ad message
			QgsMessageLog.logMessage("Started Polygon Division", MESSAGE_CATEGORY, Qgis.Info)

			# get data out of object
			# TODO: reference these properly
			layer = self.layer
			outFilePath = self.outFilePath
			target_area = self.target_area
			absorb_flag = self.absorb_flag
			direction = self.direction

			# validation that the file is projected
			if layer.crs().isGeographic():
				QgsMessageLog.logMessage("Whoops! The Polygon Divider requires a projected dataset - please save a copy with a projected CRS and try again.", MESSAGE_CATEGORY, Qgis.Critical)
				self.exception = Exception("Whoops! The Polygon Divider requires a projected dataset - please save a copy with a projected CRS and try again.")
				return False

			# used to control progress bar (only send signal for an increase)
			currProgress = 0
			self.setProgress(currProgress)

			# initial settings
			t = self.tolerance   # initial tolerance for function rooting - this is flexible now it has been divorced from the buffer
			buffer = 1e-6        # this is the buffer to ensure that an intersection occurs

			# set the direction (horizontal or vertical)
			if direction < 2:
				horizontal_flag = True
			else:
				horizontal_flag = False

			# True = cut from bottom / left, False = cut from top / right
			if (direction == 0 or direction == 2):
				forward_flag = True
			else:
				forward_flag = False

			# this is used to make sure we don't hit an insurmountable error and just repeatedly change direction
			# effectively, there are 4 directions in which we can cut, if all fail then all of these will be true and we give up
			ERROR_FLAG_0 = False	# tracks if increasing number of subdivisions failed
			ERROR_FLAG_1 = False	# tracks decreasing number of subdivisions to try and work around an error
			ERROR_FLAG_2 = False	# tracks change direction from forward to backward (or vice versa) by switching forward_flag
			ERROR_FLAG_3 = False	# tracks change cutline from horizontal to backward (or vice versa) by switching horizontal_flag

			# get fields from the input shapefile
			fieldList = self.layer.fields()

			# TODO: NEED TO CHECK IF THEY ALREADY EXIST
			# add new fields for this tool
			fieldList.append(QgsField('POLY_ID', QVariant.Int))
			fieldList.append(QgsField('UNIQUE_ID', QVariant.String))
			fieldList.append(QgsField('AREA', QVariant.Double))
			fieldList.append(QgsField('POINTX', QVariant.Int))
			fieldList.append(QgsField('POINTY', QVariant.Int))

			# create a new shapefile to write the results to
			transform_context = QgsProject.instance().transformContext()
			save_options = QgsVectorFileWriter.SaveVectorOptions()
			save_options.driverName = "ESRI Shapefile"
			save_options.fileEncoding = "UTF-8"
			writer = QgsVectorFileWriter.create(outFilePath, fieldList, QgsWkbTypes.Polygon, layer.crs(), transform_context, save_options)
			if writer.hasError() != QgsVectorFileWriter.NoError:
				QgsMessageLog.logMessage(f"Error when creating shapefile: {writer.errorMessage()}", MESSAGE_CATEGORY, Qgis.Critical)


			# define this to ensure that it's global
			subfeatures = []

			# init feature counter (for ID's)
			j = 0

			# how many sections will we have (for progress bar)
			iter = layer.getFeatures()
			totalArea = 0
			for feat in iter:
				totalArea += feat.geometry().area()
			totalDivisions = totalArea // target_area

			# check if you've been killed
			if self.isCanceled():
				# raise UserAbortedNotification('USER Killed')
				return False

			# loop through all of the features in the input data
			iter = layer.getFeatures()
			for feat in iter:

				# verify that it is a polygon
				if feat.geometry().wkbType() in [QgsWkbTypes.Polygon, QgsWkbTypes.PolygonZ,
					QgsWkbTypes.PolygonM, QgsWkbTypes.PolygonZM, QgsWkbTypes.MultiPolygon,
					QgsWkbTypes.MultiPolygonZ, QgsWkbTypes.MultiPolygonM, QgsWkbTypes.MultiPolygonZM]:

					# get the attributes to write out
					currAttributes = feat.attributes()

					# extract the geometry and sort out self intersections etc. with a buffer of 0m
					bufferedPolygon = feat.geometry().buffer(0, 15)

					# if the buffer came back as None, skip
					if bufferedPolygon is None:
						QgsMessageLog.logMessage("A polygon could not be buffered by QGIS, ignoring", MESSAGE_CATEGORY, Qgis.Info)
						continue

					# make multipolygon into list of polygons...
					subfeatures = []
					if bufferedPolygon.isMultipart():
						multiGeom = QgsGeometry()
						multiGeom = bufferedPolygon.asMultiPolygon()
						for i in multiGeom:
							subfeatures.append(QgsGeometry().fromPolygonXY(i))
					else:
						# ...OR load the feature into a list of one (it may be extended in the course of splitting if we create noncontiguous offcuts) and loop through it
						subfeatures.append(bufferedPolygon)

					#loop through the geometries
					for poly in subfeatures:

						# how many polygons are we going to have to chop off?
						nPolygons = int(poly.area() // target_area)

						# (needs to be at least 1...)
						if nPolygons == 0:
							nPolygons = 1

						# adjust the targetArea to reflect absorption if required
						if absorb_flag:
							targetArea = target_area + ((poly.area() % target_area) / nPolygons)
						else:
							targetArea = target_area

						# work out the size of a square with area = targetArea if required
						sq = targetArea**0.5

						# until there is no more dividing to do...
						while poly.area() > targetArea + t:

							# the bounds are used for the interval
							boundsR = poly.boundingBox()
							bounds = [boundsR.xMinimum(), boundsR.yMinimum(), boundsR.xMaximum(), boundsR.yMaximum()]

							# get interval and fixed coordinates (buffer otherwise there won't be an intersection between polygon and cutline!)
							if horizontal_flag:
								interval = bounds[1] + buffer, bounds[3] - buffer
								fixedCoords = bounds[0], bounds[2]
							else:
								interval = bounds[0] + buffer, bounds[2] - buffer
								fixedCoords = bounds[1], bounds[3]

							# is the interval larger than the required square? (We know the required area is > target+t because of the while statement)
							if (interval[1]-interval[0]) > sq:

								# this is the resulting area of a slice the width of sq from the polygon
								if forward_flag:
									sqArea = self.getSliceArea(interval[0] + sq - buffer, poly, fixedCoords[0], fixedCoords[1], horizontal_flag, forward_flag)	# cutting from bottom/left
								else:
									sqArea = self.getSliceArea(interval[1] - sq + buffer, poly, fixedCoords[0], fixedCoords[1], horizontal_flag, forward_flag)	# cutting from top/right

								# what is the nearest number of subdivisions of targetArea that could be extracted from that slice?
								nSubdivisions = int(round(sqArea / targetArea))

								# if the answer is 0, make it 1...
								if nSubdivisions == 0:
									nSubdivisions = 1

								# make a backup copy to reset if we move from ERROR_FLAG_0 to ERROR_FLAG_1
								nSubdivisions2 = nSubdivisions

								'''now use Brent's method to find the optimal coordinate in the variable dimension (e.g. the y coord for a horizontal cut)'''

								# if it fails, try increasing nSubdivisions (k) until it works or you get a different error
								while True:

									# how big must the target area be to support this many subdivisions?
									initialTargetArea = nSubdivisions * targetArea

									# try to split using this new value

									# try to zero the equation
									result = self.brent(interval[0], interval[1], 1e-6, t, 500, poly, fixedCoords[0], fixedCoords[1], initialTargetArea, horizontal_flag, forward_flag)

									# if it worked (no exception raised) then exit this while loop and carry on
									if result:
										break
									# otherwise, there must be an exception
									else:
										# is it a W condition error?
										if self.exception.value == "Bracket is smaller than tolerance.":

											# ...increase number of subdivisions and go around again
											nSubdivisions += 1
											continue

										# not a W condition error
										else:
											# set flag and stop trying to adjust nSubdivisions
											ERROR_FLAG_0 = True
											break

								# if that didn't work, try decreasing instead of increasing
								if ERROR_FLAG_0:

									# log message
									QgsMessageLog.logMessage("Increasing number of subdivisions didn't work, try decreasing... (Division)", MESSAGE_CATEGORY, Qgis.Warning)

									nSubdivisions = nSubdivisions2	# reset
									limit = 1
									while nSubdivisions >= limit:

										# set the flag if it's the last time around
										if nSubdivisions == limit:
											ERROR_FLAG_1 = True

										# how big must the target area be to support this many subdivisions?
										initialTargetArea = nSubdivisions * targetArea

										# try to split using this new value

										# try to zero the equation
										result = self.brent(interval[0], interval[1], 1e-6, t, 500, poly, fixedCoords[0], fixedCoords[1], initialTargetArea, horizontal_flag, forward_flag)

										# if it worked (no exception raised) then exit this while loop and carry on
										if result:
											break
										# otherwise, there must be an exception
										else:
											# ...increase number of subdivisions and go around again
											nSubdivisions -= 1
											continue

								# if increasing the subdivision size didn't help, then start trying shifting directions
								if ERROR_FLAG_1:

									# these need resetting here otherwise it won't try to cut again, just skip to the next error!
									ERROR_FLAG_0 = False
									ERROR_FLAG_1 = False

									# log message
									QgsMessageLog.logMessage("Decreasing number of subdivisions didn't work, try playing with direction... (Division)", MESSAGE_CATEGORY, Qgis.Warning)

									# switch the movement direction
									if ERROR_FLAG_2 == False:

										# log that this has been tried
										ERROR_FLAG_2 = True
										QgsMessageLog.logMessage("Reversing movement direction (Division)", MESSAGE_CATEGORY, Qgis.Warning)

										# reverse the direction of movement and try again
										forward_flag = not forward_flag
										continue

									# if the above didn't work, switch the direction of the cutline
									elif ERROR_FLAG_3 == False:

										# un-log 2, meaning that it will run again and so try the 4th direction
										ERROR_FLAG_2 = False

										# log that this has been tried
										ERROR_FLAG_3 = True
										QgsMessageLog.logMessage("Reversing cutline direction (Division)", MESSAGE_CATEGORY, Qgis.Warning)

										# reverse the cutline direction and try again
										horizontal_flag = not horizontal_flag
										continue

									# if none of the above worked, just skip it and move to a new feature
									else:
										''' WRITE THE UNSPLITTABLE POLYGON TO THE SHAPEFILE ANYWAY '''

										# make a feature with the right schema
										fet = QgsFeature()
										fet.setFields(fieldList)

										# populate inherited attributes
										for a in range(len(currAttributes)):
											fet[a] = currAttributes[a]

										# populate new attributes
										fet.setAttribute('POLY_ID', j)
										fet.setAttribute('UNIQUE_ID', str(uuid4()))
										fet.setAttribute('AREA', poly.area())

										# add the geometry to the feature
										fet.setGeometry(poly)

										# write the feature to the out file
										writer.addFeature(fet)

										# increment feature counter and
										j+=1

										# update progress bar if required
										if j // totalDivisions * 100 > currProgress:
											self.setProgress(j // totalDivisions * 100)

										# log that there was a problem
										QgsMessageLog.logMessage("There was an un-dividable polygon in this dataset.", MESSAGE_CATEGORY, Qgis.Warning)

										# on to the next one, hopefully with more luck!
										continue

								# if it worked, reset the flags
								ERROR_FLAG_0 = False
								ERROR_FLAG_1 = False
								ERROR_FLAG_2 = False
								ERROR_FLAG_3 = False

								# create the desired cutline as lists of QgsPoints
								if horizontal_flag:
									line = [(fixedCoords[0], result), (fixedCoords[1], result)] # horizontal split
								else:
									line = [(result, fixedCoords[0]), (result, fixedCoords[1])] # vertical split

								# calculate the resulting polygons - poly will be sliced again, initialSlice will be subdivided
								poly, initialSlice, residuals = self.splitPoly(poly, line, horizontal_flag, forward_flag)

								# put the residuals in the list to be processed
								subfeatures += residuals

							# bounds not bigger than sq, so no division necessary, just subdivide this last one directly (nothing will happen if it can't be subdivided)
							else:
								# set the remainder of the polygon as the final slice, and poly to an empty polygon
								initialSlice = poly
								poly = QgsGeometry.fromPolygonXY([[]])

								# what is the nearest number of subdivisions of targetArea that could be extracted from that slice? (must be at least 1)
								# TODO: verify this doesn't need rounding
								nSubdivisions = int(initialSlice.area() // targetArea) # shouldn't need rounding...
								if nSubdivisions == 0:
									nSubdivisions = 1

							#...then divide that into sections of targetArea
							for k in range(nSubdivisions-1):	# nCuts = nPieces - 1

								# the bounds are used for the interval
								sliceBoundsR = initialSlice.boundingBox()
								sliceBounds = [sliceBoundsR.xMinimum(), sliceBoundsR.yMinimum(), sliceBoundsR.xMaximum(), sliceBoundsR.yMaximum()]

								# get the slice direction (opposite to main direction)
								sliceHorizontal = not horizontal_flag

								if sliceHorizontal:
									# get interval and fixed coordinates
									sliceInterval = sliceBounds[1] + buffer, sliceBounds[3] - buffer # buffer otherwise there won't be an intersection between polygon and cutline!
									sliceFixedCoords = sliceBounds[0], sliceBounds[2]
								else:
									# get interval and fixed coordinates
									sliceInterval = sliceBounds[0] + buffer, sliceBounds[2] - buffer # buffer otherwise there won't be an intersection between polygon and cutline!
									sliceFixedCoords = sliceBounds[1], sliceBounds[3]

								# restore the tolerance (may be adjusted in the below loop)
								tol = t

								# infinite loop
								while True:

									# brent's method to find the optimal coordinate in the variable dimension (e.g. the y coord for a horizontal cut)

									# search for result
									sliceResult = self.brent(sliceInterval[0], sliceInterval[1], 1e-6, tol, 500, initialSlice, sliceFixedCoords[0], sliceFixedCoords[1], targetArea, sliceHorizontal, forward_flag)

									# stop searching if result is found
									if sliceResult:
										break
									else:
										# if it is a W condition error, double the tolerance
										if self.exception.value == "Bracket is smaller than tolerance.":
											QgsMessageLog.logMessage(self.exception.value + ": increasing tolerance (Subdivision)", MESSAGE_CATEGORY, Qgis.Warning)

											# double the tolerance and try again
											tol *= 2
											continue

										# otherwise, give up and try something else
										else:
											# set the flag that this has been tried and move on
											ERROR_FLAG_1 = True
											break

								''' IF THE ABOVE DIDNT WORK THEN WE NEED TO TRY MORE DRASTIC MEASURES '''

								# try reversing the movement direction
								if ERROR_FLAG_1 and not ERROR_FLAG_2: # (NB: Subdivision does not use Errorflag 0)

									# log that this has been tried
									ERROR_FLAG_2 = True
									QgsMessageLog.logMessage("Reversing movement direction (Subdivision)", MESSAGE_CATEGORY, Qgis.Warning)

									# reverse the direction of movement and try again
									forward_flag = not forward_flag
									continue

								# if that didn't work, switch the direction of the cutline
								elif ERROR_FLAG_1 and not ERROR_FLAG_3:

									# log that this has been tried
									ERROR_FLAG_3 = True
									QgsMessageLog.logMessage("Reversing cutline direction (Subdivision):", MESSAGE_CATEGORY, Qgis.Warning)

									# reverse the cutline direction and pass back to the outer division to try again in the opposite direction (otherwise we would get long thin strips, not squares)
									horizontal_flag = not horizontal_flag
									break	# this should mean that the 'else' for this statement will never be reached

								# if it worked, reset the flags
								ERROR_FLAG_1 = False
								ERROR_FLAG_2 = False
								ERROR_FLAG_3 = False

								# create the desired cutline as lists of QgsPoints
								if horizontal_flag:
									sliceLine = [(sliceResult, sliceFixedCoords[0]), (sliceResult, sliceFixedCoords[1])]	# horizontal split
								else:
									sliceLine = [(sliceFixedCoords[0], sliceResult), (sliceFixedCoords[1], sliceResult)]	# vertical split

								# calculate the resulting polygons - initialSlice becomes left (to be chopped again)
								initialSlice, right, residuals = self.splitPoly(initialSlice, sliceLine, sliceHorizontal, forward_flag)

								# put the residuals in the list to be processed
								subfeatures += residuals

								''' WRITE TO SHAPEFILE '''

								# make a feature with the right schema
								fet = QgsFeature()
								fet.setFields(fieldList)

								# populate inherited attributes
								for a in range(len(currAttributes)):
									fet[a] = currAttributes[a]

								# calculate representative point
								pt = right.pointOnSurface().asPoint()

								# populate new attributes
								fet.setAttribute('POLY_ID', j)
								fet.setAttribute('UNIQUE_ID', str(uuid4()))
								fet.setAttribute('AREA', right.area())
								fet.setAttribute('POINTX', pt[0])
								fet.setAttribute('POINTY', pt[1])

								# add the geometry to the feature
								fet.setGeometry(right)

								# write the feature to the out file
								writer.addFeature(fet)

								# increment feature counter and
								j+=1

								# update progress bar if required
								if int((j*1.0) / totalDivisions * 100) > currProgress:
									currProgress = int((j*1.0) / totalDivisions * 100)
									self.setProgress(currProgress)

							## WRITE ANY OFFCUT FROM SUBDIVISION TO SHAPEFILE

							# make a feature with the right schema
							fet = QgsFeature()
							fet.setFields(fieldList)

							# populate inherited attributes
							for a in range(len(currAttributes)):
								fet[a] = currAttributes[a]

							# calculate representative point
							pt = initialSlice.pointOnSurface().asPoint()

							# populate new attributes
							fet.setAttribute('POLY_ID', j)
							fet.setAttribute('UNIQUE_ID', str(uuid4()))
							fet.setAttribute('AREA', initialSlice.area())
							fet.setAttribute('POINTX', pt[0])
							fet.setAttribute('POINTY', pt[1])

							# add the geometry to the feature
							fet.setGeometry(initialSlice)

							# write the feature to the out file
							writer.addFeature(fet)

							# increment feature counter and
							j+=1

							# update progress bar if required
							if int((j*1.0) / totalDivisions * 100) > currProgress:
								currProgress = int((j*1.0) / totalDivisions * 100)
								self.setProgress(currProgress)

						try:

							## WRITE  ANY OFFCUT FROM DIVISION TO SHAPEFILE

							# make a feature with the right schema
							fet = QgsFeature()
							fet.setFields(fieldList)

							# populate inherited attributes
							for a in range(len(currAttributes)):
								fet[a] = currAttributes[a]

							# calculate representative point
							pt = poly.pointOnSurface().asPoint()

							# populate new attributes
							fet.setAttribute('POLY_ID', j)
							fet.setAttribute('UNIQUE_ID', str(uuid4()))
							fet.setAttribute('AREA', poly.area())
							fet.setAttribute('POINTX', pt[0])
							fet.setAttribute('POINTY', pt[1])

							# add the geometry to the feature
							fet.setGeometry(poly)

							# write the feature to the out file
							writer.addFeature(fet)

							# increment feature counter and
							j+=1

							# update progress bar if required
							if int((j*1.0) / totalDivisions * 100) > currProgress:
								currProgress = int((j*1.0) / totalDivisions * 100)
								self.setProgress(currProgress)

						except:
							# this just means that there is no offcut, which is no problem!
							pass
				else:
					QgsMessageLog.logMessage("Whoops! That dataset isn't polygons!", MESSAGE_CATEGORY, Qgis.Critical)
					# raise Exception("Whoops! That dataset isn't polygons!")
					self.exception = Exception("Whoops! That dataset isn't polygons!")
					return False

			if self.isCanceled():
				# raise UserAbortedNotification('USER Killed')
				return False

			# finally, open the resulting file and return it
			layer = QgsVectorLayer(outFilePath, 'Divided Polygon', 'ogr')
			if layer.isValid():
				# return success
				return True
			else:
				self.exception = Exception("Output Dataset Invalid")
				return False

		# catch-all to prevent exceptions being buried by QGIS
		except Exception as e:
			self.exception = e
			return False

class PolygonDivider:
	"""QGIS Plugin Implementation."""

	def __init__(self, iface):
		"""Constructor.

		:param iface: An interface instance that will be passed to this class
			which provides the hook by which you can manipulate the QGIS
			application at run time.
		:type iface: QgsInterface
		"""
		# Save reference to the QGIS interface
		self.iface = iface
		# initialize plugin directory
		self.plugin_dir = os.path.dirname(__file__)
		# initialize locale
		locale = QSettings().value('locale/userLocale')[0:2]
		locale_path = os.path.join(
			self.plugin_dir,
			'i18n',
			'PolygonDivider_{}.qm'.format(locale))

		if os.path.exists(locale_path):
			self.translator = QTranslator()
			self.translator.load(locale_path)
			QCoreApplication.installTranslator(self.translator)

		# Declare instance attributes
		self.actions = []
		self.menu = self.tr(u'&Polygon Divider')

		# JJH: declare task manager to handle the backcround task
		self.tm = QgsApplication.taskManager()

		# Check if plugin was started the first time in current QGIS session
		# Must be set in initGui() to survive plugin reloads
		self.first_start = None


	# noinspection PyMethodMayBeStatic
	def tr(self, message):
		"""Get the translation for a string using Qt translation API.

		We implement this ourselves since we do not inherit QObject.

		:param message: String for translation.
		:type message: str, QString

		:returns: Translated version of message.
		:rtype: QString
		"""
		# noinspection PyTypeChecker,PyArgumentList,PyCallByClass
		return QCoreApplication.translate('PolygonDivider', message)


	def add_action(
		self,
		icon_path,
		text,
		callback,
		enabled_flag=True,
		add_to_menu=True,
		add_to_toolbar=True,
		status_tip=None,
		whats_this=None,
		parent=None):
		"""Add a toolbar icon to the toolbar.

		:param icon_path: Path to the icon for this action. Can be a resource
			path (e.g. ':/plugins/foo/bar.png') or a normal file system path.
		:type icon_path: str

		:param text: Text that should be shown in menu items for this action.
		:type text: str

		:param callback: Function to be called when the action is triggered.
		:type callback: function

		:param enabled_flag: A flag indicating if the action should be enabled
			by default. Defaults to True.
		:type enabled_flag: bool

		:param add_to_menu: Flag indicating whether the action should also
			be added to the menu. Defaults to True.
		:type add_to_menu: bool

		:param add_to_toolbar: Flag indicating whether the action should also
			be added to the toolbar. Defaults to True.
		:type add_to_toolbar: bool

		:param status_tip: Optional text to show in a popup when mouse pointer
			hovers over the action.
		:type status_tip: str

		:param parent: Parent widget for the new action. Defaults None.
		:type parent: QWidget

		:param whats_this: Optional text to show in the status bar when the
			mouse pointer hovers over the action.

		:returns: The action that was created. Note that the action is also
			added to self.actions list.
		:rtype: QAction
		"""

		# Create the dialog (after translation) and keep reference
		self.dlg = PolygonDividerDialog()

		icon = QIcon(icon_path)
		action = QAction(icon, text, parent)
		action.triggered.connect(callback)
		action.setEnabled(enabled_flag)

		if status_tip is not None:
			action.setStatusTip(status_tip)

		if whats_this is not None:
			action.setWhatsThis(whats_this)

		if add_to_toolbar:
			# Adds plugin icon to Plugins toolbar
			self.iface.addToolBarIcon(action)

		if add_to_menu:
			self.iface.addPluginToVectorMenu(
				self.menu,
				action)

		self.actions.append(action)

		return action


	def select_output_file(self):
		"""
		* Open file browser
		"""

		# get filename from dialog
		filename = QFileDialog.getSaveFileName(self.dlg, "Select output file ","", '*.shp')[0]

		# verify that a name was selected
		if filename != "":

			# clear previous value
			self.dlg.lineEdit_2.clear()

			# make sure that an extension was included
			if filename[-4:] != '.shp':
				filename += '.shp'

			# put the result in the text box on the dialog
			self.dlg.lineEdit_2.setText(filename)


	def radio_handler(self):
		"""
		* Respond to radio button toggle by setting which text box is active and
		*  setting the flag
		"""
		# see which is checked and enable / disable the text boxes accordingly
		if self.dlg.radioButton_2.isChecked():
			self.dlg.lineEdit.setDisabled(True)
			self.dlg.lineEdit_3.setEnabled(True)
		elif self.dlg.radioButton.isChecked():
			self.dlg.lineEdit.setEnabled(True)
			self.dlg.lineEdit_3.setDisabled(True)


	def initGui(self):
		"""Create the menu entries and toolbar icons inside the QGIS GUI."""

		icon_path = ':/plugins/polygon_divider/icon.png'
		self.add_action(
			icon_path,
			text=self.tr(u'Polygon Divider'),
			callback=self.run,
			parent=self.iface.mainWindow())

		# will be set False in run()
		self.first_start = True


	def unload(self):
		"""Removes the plugin menu item and icon from QGIS GUI."""
		for action in self.actions:
			self.iface.removePluginVectorMenu(
				self.tr(u'&Polygon Divider'),
				action)
			self.iface.removeToolBarIcon(action)


	def run(self):
		"""
		* Run method that performs all the real work
		"""

		# Create the dialog with elements (after translation) and keep reference
		# Only create GUI ONCE in callback, so that it will only load when the plugin is started
		if self.first_start == True:

			# toggle flag
			self.first_start = False

			# get reference to the dialog object
			self.dlg = PolygonDividerDialog()

			''' THIS SETUP FIRST TIME IT LOADS ONLY '''

			# button listener - launch file browser for output file
			self.dlg.pushButton.clicked.connect(self.select_output_file)

			# radio button listeners - flip enabled field
			self.dlg.radioButton.toggled.connect(self.radio_handler)
			self.dlg.radioButton_2.toggled.connect(self.radio_handler)

		''' THIS SETUP EVERY TIME IT LOADS '''

		# populate comboBox with the active layers
		self.dlg.comboBox.clear()	# need to clear here or it will add them all again every time the dialog is opened
		layers = QgsProject.instance().layerTreeRoot().findLayers()
		self.dlg.comboBox.addItems([layer.name() for layer in layers])

		# populate comboBox_2 with the possible directions
		self.dlg.comboBox_2.clear() # need to clear here or it will add them all again every time the dialog is opened
		self.dlg.comboBox_2.addItems(['left to right', 'right to left', 'bottom to top', 'top to bottom'])

		# show the dialog
		self.dlg.show()

		# Run the dialog event loop
		result = self.dlg.exec_()

		# See if OK was pressed
		if result:

			''' RUN THE TOOL '''

			# get user settings
			inFile = layers[self.dlg.comboBox.currentIndex()].layer()
			outFilePath = self.dlg.lineEdit_2.text()
			targetArea = float(self.dlg.lineEdit.text())
			absorbFlag = self.dlg.checkBox.isChecked()
			direction = self.dlg.comboBox_2.currentIndex()
			tolerance = float(self.dlg.lineEdit_4.text())

			# if the user has selected number of features option
			if self.dlg.radioButton_2.isChecked():

				# sum the area of all features in the input layer
				total_area = 0
				for feature in inFile.getFeatures():
					total_area += feature.geometry().area()

				# calculate the desired target area
				targetArea = int((total_area / float(self.dlg.lineEdit_3.text())) + 0.5)

			# launch the task
			self.tm.addTask(PolygonDividerTask(inFile, outFilePath, targetArea, absorbFlag, direction, tolerance))