Quarter Geometries

quart_geom.py

import arcpy
import os
import io
import sys


def quart_geom(geom):
    """
    Quarter an arcpy geometry
    :param geom: arcpy.geometry
    :return: list
    """

    # Get the extent of the input geometry, get corner points and calculte X/Y coordinates
    min_x = geom.extent.XMin
    max_x = geom.extent.XMax
    mean_x = (max_x + min_x) / 2
    min_y = geom.extent.YMin
    max_y = geom.extent.YMax
    mean_y = (max_y + min_y) / 2

    # construct list of new corner points for the four quarters
    quart1_coords = [[min_x, min_y], [mean_x, min_y], [mean_x, mean_y], [min_x, mean_y], [min_x, min_y]]
    quart2_coords = [[mean_x, min_y], [max_x, min_y], [max_x, mean_y], [mean_x, mean_y], [mean_x, min_y]]
    quart3_coords = [[mean_x, mean_y], [max_x, mean_y], [max_x, max_y], [mean_x, max_y], [mean_x, mean_y]]
    quart4_coords = [[min_x, mean_y], [mean_x, mean_y], [mean_x, max_y], [min_x, max_y], [min_x, mean_y]]

    quarts_coords = [quart1_coords, quart2_coords, quart3_coords, quart4_coords]

    # construct polygon geometries, using quart coordinate lists and add them to a list
    quarts = list()
    for quart_coords in quarts_coords:
        # Create a Polygon object based on the array of points
        # Append to the list of Polygon objects
        quarts.append(
            arcpy.Polygon(
                arcpy.Array([arcpy.Point(*coords) for coords in quart_coords])))

    # clip input feature using newly created quarts and add them to a list
    geom_quarts = list()
    for quart in quarts:
        try:
            #if geom.overlaps(quart):
                geom_quarts.append(geom.clip(quart.extent))
        except ValueError:
            # skip if something goes wrong - not ideal - have to find a better way to handle this
            # this happens a lot with very large geometries. Seems to be an Arcpy issue?
            print "skip"
            
    # return list containing up to four geometries
    return geom_quarts

def quart_geom_max_byte(geometry, max_bytes=1e+7):
    """
    This function quarters an arcpy geometry.
    If WKT of resulting quarter is larger than given max_bytes it will be quartered again and again
    until threshold is met.
    Return list of quartered arcpy geometries
    :param geometry: arcpy.geometry
    :param max_bytes: float
    :return: list
    """

    passed_geoms = list()
    large_geoms = list()
    done = False

    while not done:
        done = True

        # For the first run add input geometry to check_geoms as a list
        if not passed_geoms and not large_geoms:
            check_geoms = [geometry]
        # otherwise use the output results from the geometries quartered in the previous run
        # and reset large_geoms
        else:
            check_geoms = large_geoms
            large_geoms = list()

        # fall all geometries not yet checked, test if WKT bytes size is larger than threshold
        # if yes, quarter it, if not append it to passed_geoms
        for geom in check_geoms:
            try:
                if sys.getsizeof(geom.WKT) > max_bytes:
                    large_geoms += quart_geom(geom)
                    done = False
                else:
                    passed_geoms.append(geom)
            except SystemError:
                # geom.WKT fails for very large geometries
                # in that case, just assume that it is to big and quarter it

                print "too large"
                large_geoms += quart_geom(geom)
                done = False

    # return list of quartered geoms
    return passed_geoms


if __name__ == "__main__":

    # Get GDAM feature class and select relevant fields
    fc = r"D:\GIS Data\Global\gadm28.gdb\gadm"

    cur = arcpy.da.SearchCursor(fc, ["ISO", "NAME_0", "HASC_1", "NAME_1",
                                     "HASC_2", "NAME_2", "HASC_3", "NAME_3",
                                     "NAME_4", "NAME_5", "SHAPE@"],
                                sql_clause=(None, 'ORDER BY ISO'))

    # Define output path (needs to be empty!!)
    gadm_path = r"D:\temp\gadm\tsv"

    # Loop through all rows and convert write them to TSV file
    # if geometry is too large, quarter it
    counter = 0
    for row in cur:
        # quarter geometries
        geoms = quart_geom_max_byte(row[10], max_bytes=1e+7)
        if len(geoms) > 1:
            print len(geoms)
            #print u" ".join(i.encode('utf8') if isinstance(i, basestring) else u"" for i in row[:-1])

        # define output file name using country code
        # one file for each country
        gadm_file = os.path.join(gadm_path, "{}.tsv".format(row[0]))
        
        #if output file does not yet exist we are dealing with a new country
        if not os.path.exists(gadm_file):
            # if this is not the first loop
            # close previous file
            if counter > 0:
                gadm.close()
                counter = 0
            # Create new file
            gadm = io.open(gadm_file, 'w', encoding='utf-8')
            print gadm_file
            #Had some encoding issue. This might now actually work without encode/decode. Seems now redundant but I haven't tested it without...

        # write row to currently opened file
        # if geometry was quartered several rows with the same attributs but different geometries will be added
        #tsv_rows = list()
        for geom in geoms:
            if counter > 0:
                gadm.write(u"\n")
            counter += 1
            if geom.WKT != "MULTIPOLYGON EMPTY":
                new_row = row[:-1] + (geom.WKT,)
                #tsv_rows.append(u"\t".join(i.encode('utf8').decode('UTF-8') if isinstance(i, basestring) else u"" for i in new_row)])
                gadm.write(u"\t".join(i if isinstance(i, basestring) else u"" for i in new_row))
        #gadm.write(u"\n".join(tsv_rows))

    # Finally, close last file
    gadm.close()