After buildings have been extracted or manually digitized from satellite imageries, they often produce building outlines that are not correctly orthogonal at their corners whereas they are in reality. It is then important to regularize these outlines for further use for mapmaking, further analysis and eventually, decision making.

The goal of the project is to create new buildings footprint datasets that enforce rectangular corners. That is, wherever the original corner angles fall between 90°- ε and 90°+ ε. Where ε is an inaccuracy measure.

What you will need?

• Basic knowledge of Python (Vanilla & Numpy) and GIS.
• ArcPy (ArcGIS 10.7) -- Python 2 is the supported version
• Choice code editor (I use VSCode 😉)
• Sample dataset (projected in UTM) has been provided in a Github repository.

For simplicity, the context of valid polygons is used as input. A valid polygon, in this case, has unambiguous exteriors, no touching segments except at vertices, non-zero length, the clockwise direction of outer rings, non-zero area and with no overlap. That is, topologically simple shapes with no multi-parts.

The geometry below is an example context considered.

Workflow

The workflow developed has been modularized using functional programming approaches. Before we go further, it is important to highlight that input datasets are best in UTM projection as Earth distance and angle measurements will be needed.

Utility functions

The function below validates the file extension.

def controlExtension(inName, ext):
    if inName.rfind('.') > 0:
        return inName[:inName.find('.')] + ext
    else:
        return inName + ext

Now, we check to see if the input data exists.

def checkExistence(pathList):
    check = True
    for data in pathList:
        if not arcpy.Exists(data):
            check = False
            print '! dataset ' + data + ' is missing'
            break
    return check

Creating subdirectories in workspace to store temporary files as well as outputs.

def createSubdir(workspace, subdirList):
    for subdir in subdirList:
        if not os.path.isdir(workspace + '/' + subdir):
            os.mkdir(os.path.join(workspace, subdir))

We also need a function that creates a complete path for processing input, temporary and output files.

def completePath(workspace, subdir, nameList):
    for ix in range(len(nameList)):
        nameList[ix] = workspace + '/' + subdir + '/' + nameList[ix]
    return nameList

And we also need to specify the file extension of the output datasets which will eventually be used for mapmaking and decision-making.

def outputFiles(outputList):
    file = []
    for output in outputList:
        file.append(str(output) + '.shp')
    return file

From time to time, we will be storing and re-using values from our computation in the attribute table fields (i.e. columns) and so we will need a utility function to retrieve the field names.

def getFieldNames(table):
    fnames = []
    fields = arcpy.ListFields(table)
    if fields:
        for field in fields:
            fnames.append(field.name)
    return fnames

Pre-processing functions

One common situation in manually digitized features is redundant vertices. These vertices are regarded as redundant because they could have been just a single line segment with a start and end vertex. The below function reuses the SimplifyPolygon_cartography. The algorithm specified in the "POINT_REMOVE" which is based on the Douglas-Peucker Algorithm in simplifying geometries.

def simplifyBuilding(inFC, outFC):
    arcpy.SimplifyPolygon_cartography(inFC, outFC, algorithm="POINT_REMOVE", tolerance=1, minimum_area=0)

Similarly, simplifying the features also means that we should have just the needed outlines as sides of a building corner. Therefore, we then proceed to break the polygon into lines using a pre-defined function in ArcPy called PolygonToLine_management.

def polygonToLine(inFC, outFC):
    arcpy.PolygonToLine_management(inFC, outFC)

We proceed to split the lines from the last step at their vertices using yet another available function, SplitLine_management from ArcPy.

def polylineToSegments(inFC, outFC):
    arcpy.SplitLine_management(inFC, outFC)

Processing functions

There are several ways to enforce that the four corners of a building are at 90° as they are in reality. The approach selected is to adjust the three sides of a building polygon relative to a stable dominant orientation of an exterior wall. The dominant segment would be the longest side.

To begin, we get each building coordinates in arrays using the following function.

def _building_arr_(geometry):
    """Return building coordinates."""
    def _split_(part):
        yield [(p.X, p.Y) for p in part if p]

    coords =[]
        output = []
        w = np.where(np.in1d(part, None, invert=False))[0]
        bits = np.split(part, w)
        for bit in bits:
            geom_sub = _split_(bit)
            output.append(np.array(*geom_sub).squeeze())
        coords.append(np.asarray(output).squeeze())
    return np.asarray(coords)

We also need to get the angle at each vertex, that is, where two line segments meet. The angles are stored in the attribute table as type "Double".

def buildingVertexAngle(table, angles):

    arcpy.AddField_management(table,"Angle","DOUBLE")

    pointer = 0

    with arcpy.da.UpdateCursor(table, 'Angle') as uCur:
        orient = []
        for angle in angles:
            for i in angle:
                orient.append(i)

        for row in uCur:
            row[0] = orient[pointer]

            pointer += 1
            uCur.updateRow(row)

        return orient

Likewise, we store the error, ε , 90°+/- ε in another field of our attributes table.

def buildingVertexError(table, angles):
    arcpy.AddField_management(table,"Angle_Err","DOUBLE")

    pointer = 0

    with arcpy.da.UpdateCursor(table, 'Angle_Err') as uCur:
        angleError = []
        for angle in angles:
            error = 90 - angle
            angleError.append( error )

        for row in uCur:
            row[0] = angleError[pointer]

            pointer += 1
            uCur.updateRow(row)

        return angleError

We then proceed to derive the length of the building exterior walls and add same as a field to the attribute table for later use.

def buildingLengths(table, lengths):

    arcpy.AddField_management(table,"Length","DOUBLE")

    pointer = 0

    with arcpy.da.UpdateCursor(table, 'Length') as uCur:
        buildLengths = []                                  
        for length in lengths:
            for i in length:
                buildLengths.append(i)

        for row in uCur:
            row[0] = buildLengths[pointer]

            pointer += 1
            uCur.updateRow(row)
        return buildLengths

Because I have decided to approach enforcing the orthogonality of the building walls by making the longest side the reference side to adjust other sides, the next thing is to identify the longest side.

def geoMaxLength(inFC):
    coords = geomCoords(inFC)
    groups = [] 
    buildings = dict()
    items = 0

    for i, coord in enumerate(coords):
        groups.append(coord)

        if i % 4 == 0:
            groups = []
            groups.append(coords[i])

        if len(groups) == 4:
            lengths = []
            for j in groups:
                length = _length_(j)
                lengths.append(length)
            lengths = np.array(lengths)

            maxlen = np.max(lengths) 
            loc = np.argmax(lengths)

            buildings[items] = dict()
            buildings[items]['max_coord'] = groups[loc]
            buildings[items]['max_len'] = maxlen
            buildings[items]['max_index'] = loc

            items = items + 1

    return buildings

Now we know the length of the building walls, the longest side of the building, the angle at each vertex, the angle errors and coordinates at each vertex. We then proceed to generate new coordinates at each vertex adjusted for 90° with reference to the longest side of the building.

def calculateNewCoord(a, b, alpha, min_length):     
    aa = abs(min_length * (np.cos(alpha + np.pi/2.)) - a)
    bb = abs(min_length * (np.sin(alpha + np.pi/2.)) - b)

    return a, b, aa, bb

We subsequently use the new coordinates to generate new line segments, polylines and ultimately, building polygons that will be saved as a shapefile. The coordinate system (UTM) originally used in the input datasets is re-used for the output building datasets.

Conclusion

We see that it is possible to make building polygons exactly rectangular at their corners. To put things to test, performance-wise, a number of buildings were used and the outcome can be seen below:

| Number of buildings | Time of execution  |
|---------------------|--------------------|
| 4                   | 7.98099994659 secs |
| 50                  | 22.1856946297 secs |
| 100                 | 36.2016902341 secs |

And... that's it. We now have a working approach to enforcing the orthogonality of building segments.

Thank you for reading. :)

Too long and just interested in seeing the code with loads of helpful comments and documentation, head over to my repository.