Splitting square-based polygon into smaller square areas in python

Question

I have a polygon like this in a coordinate grid: [(0, 0), (0, 1), (0, 2), (0, 3), (1, 0), (1, 1), (1, 2), (1, 3), (2, 2), (2, 3)].

Is there any way in Python 2.7 I can divide this polygon into smaller square areas of arbitrary side longitude (passed as a parameter)? basically what I need is to divide square-shaped polygons into smaller square areas.

Any help would be quite appreciated!

Answer 1

This can be done by iterating through the list of coordinates and determining if the squares are able to be made (that all of the internal coordinates are in the list of coordinates).

The code below will iterate through coordinates in the list and check if they can be made into full squares. It assumes that the list of coordinates is sorted by x then by y (where the coordinates are (x, y) pairs)

It works by adding the new values to the grid list and iterating through it.

example grid returning 2 1x1 squares:

[(0,0),(0,1),(0,2),
 (1,0),(1,1),(1,2)]

function:

import math
import numpy

def getArray(grid):
    n = numpy.zeros((grid[-1][0]+1, grid[-1][1]+1))
    for (y,x) in grid:
        n[y, x] = 1
    return n

# Determines if the new point is within the bounds
def getBoundingSquare(newCoord, npArr):
    try:
        if npArr[int(math.floor(newCoord[0])),int(math.floor(newCoord[1]))] == 1 and \
        npArr[int(math.floor(newCoord[0])),int(math.ceil(newCoord[1]))] == 1 and \
        npArr[int(math.ceil(newCoord[0])),int(math.floor(newCoord[1]))] == 1 and \
        npArr[int(math.ceil(newCoord[0])),int(math.ceil(newCoord[1]))] == 1:
            return 1
        else:
            return 0
    except IndexError:
        return 0

# Creates the new points using the desired side length
def interpolator(grid, side_length):
    startCorner = grid[0]
    endCorner = grid[-1]
    npArr = getArray(grid)
    newGrid = []
    if side_length < 1:
        exprY = int((endCorner[0]+1)*1//side_length-1)
        exprX = int((endCorner[1]+1)*1//side_length-1)
    else:
        exprY = int((endCorner[0]+1))
        exprX = int((endCorner[1]+1))
    for y in range(startCorner[0], exprY):
        for x in range(startCorner[1], exprX):
            newCoord = (y*side_length+startCorner[0], x*side_length+startCorner[1])
            newCoord2 = (float(y+startCorner[0]), float(x+startCorner[1]))
            if getBoundingSquare(newCoord, npArr):
                newGrid.append(newCoord)
            if getBoundingSquare(newCoord2, npArr) and newCoord2 not in newGrid:
                newGrid.append(newCoord2)
    newGrid.sort()
    return newGrid

def subdivide(grid, side_length):
    grid = interpolator(grid, float(side_length))
    subSquares = []
    while len(grid) >= 4:
        sy, sx = grid[0]
        if (sy+side_length, sx+side_length) in grid:
            square = []
            for y in range(2):
                for x in range(2):
                    if (sy+y*side_length, sx+x*side_length) in grid:
                        square.append((sy+y*side_length, sx+x*side_length))
                        if not(y == 1 or x == 1):
                            grid.remove((sy+y*side_length, sx+x*side_length))

            if square not in subSquares and (len(square) == (side_length+1)**2 or len(square) == 4):
                subSquares.append(square)
            (startY, startX) = square[0]
            (endY, endX) = square[-1]
            counter = 0
            while counter < len(grid):
                item = grid[counter]
                if (item[0] < endY and item[1] < endX):
                    grid.remove(item)
                else:
                    counter += 1
        else:
            grid.pop(0)
    allowed = 0
    for item in grid:
        for square in subSquares:
            if item in square:
                allowed += 1
                continue
    if len(grid) > allowed:
        print 'Could not divide entire polygon'
    for square in subSquares:
        print square
    return subSquares

This does not return overlapping squares. This requires numpy to be installed.

Answer 2

I think this code is not very efficient, but it does the trick for me. Basically, it takes a dictionary of polygons (each polygon defined by tuples of coordinates x, y), in a way that every polygon could be represented by smaller squares of length la_radius and then it divides every such square into smaller squares of length sa_radius (as requested):

def divide_big_areas(big_areas,sa_radius):
    sa_coordinates = {}
    for b_area_key, b_area_value in big_areas.iteritems():
        sa_counter = 0
        #This If below is just to be a little more efficient by recognizing polygons fromed only by an isolated square
        if len(b_area_value) == 4:
            x = min(x for x, y in b_area_value)
            y = min(y for x, y in b_area_value)
            temp_area = [(x,y),(x,y+la_radius),(x+la_radius,y),(x+la_radius,y+la_radius)]
            num_sa = int(max(x for x, y in temp_area)*max(y for x, y in temp_area)/(sa_radius**2))
            for num_area in range(num_sa):
            #creating smaller areas (squared areas of size sa_radius)
                sa_coordinates[b_area_key+'_sa'+str(sa_counter)] = [(x,y),(x+sa_radius,y),(x+sa_radius,y+sa_radius),(x,y+sa_radius)]
                sa_counter += 1
                if x + 2*sa_radius <= max(x for x, y in temp_area) and y + 2*sa_radius <= max(y for x, y in temp_area):
                    x = x + sa_radius
                elif x + 2*sa_radius >= max(x for x, y in temp_area) and y + 2*sa_radius <= max(y for x, y in temp_area):
                    y = y + sa_radius
                    x = 0
                elif x + 2*sa_radius <= max(x for x, y in temp_area) and y + 2*sa_radius >= max(y for x, y in temp_area):
                    x = x + sa_radius
                else:
                    break
        else:
            for x1, y1 in b_area_value:
                x = x1
                y = y1
                if ((x+la_radius,y) in b_area_value and (x,y+la_radius) in b_area_value and (x+la_radius,y+la_radius) in b_area_value and x+sa_radius <= x+la_radius and y+sa_radius <= y+la_radius):
                    temp_area = [(x,y),(x,y+la_radius),(x+la_radius,y),(x+la_radius,y+la_radius)]
                    num_sa = int(la_radius**2/(sa_radius**2))
                    for num_area in range(num_sa):
                        #creating smaller areas (squared areas of size sa_radius)
                        sa_coordinates[b_area_key+'_sa'+str(sa_counter)] = [(x,y),(x+sa_radius,y),(x+sa_radius,y+sa_radius),(x,y+sa_radius)]
                        sa_counter += 1
                        if x + 2*sa_radius <= max(x for x, y in temp_area) and y + 2*sa_radius <= max(y for x, y in temp_area):
                            x = x + sa_radius
                        elif x + 2*sa_radius >= max(x for x, y in temp_area) and y + 2*sa_radius <= max(y for x, y in temp_area):
                            y = y + sa_radius
                            x = x1
                        elif x + 2*sa_radius <= max(x for x, y in temp_area) and y + 2*sa_radius >= max(y for x, y in temp_area):
                            x = x + sa_radius
                        else:
                            break      
    return (sa_coordinates)