Is centering existing polygons on a point layer possible?

Question

I'm writing a python script that automates our grid maps. Everything works fine other than the occasionally awkwardly placed grids.

I've been using Create Grid Index Feature, which puts the orgin at the bottom left of the layers extent. This is fine for large datasets as typically it reflects best fit. I even have the code componsate for when only a single grid is required (center on the points layer, make extent equal to dataframe at a scale of 9000, then get xmin and ymin of data frame). The problem is the smaller to medium sized point clusters.

I'm curious if there is a way to align a grid to best fit around the points layer. Currently, a line of points could starts at the bottom left corner of the grid and go through one sheet of my grid map and have a single point on the second page. I would like to align the grid around the points so the points are evenly distributed.

Is this something that's possible within the confines of ArcPy?

Approach Utilized

Using Richard Fairhurst's answer, I realized I was going about this wrong. Using his example I decided to get the center coordinate of my grid and points layer, then I took the difference of the two centroids and utilized that value to shift all features in said layer. I stubled upon and utilized the following function to do so; which I found at ArcPy Café.

def shift_features(in_features, x_shift=None, y_shift=None):
    """
    Shifts features by an x and/or y value. The shift values are in
    the units of the in_features coordinate system.

    Parameters:
    in_features: string
    An existing feature class or feature layer.  If using a
    feature layer with a selection, only the selected features
    will be modified.

    x_shift: float
    The distance the x coordinates will be shifted.

    y_shift: float
    The distance the y coordinates will be shifted.
    """
    with arcpy.da.UpdateCursor(in_features, ['SHAPE@XY']) as cursor:
        for row in cursor:
            cursor.updateRow([[row[0][0] + (x_shift or 0), row[0][1] + (y_shift or 0)]])
    return

http://gis.stackexchange.com/questions/154975/how-to-minimise-number-of-dynamic-pages-to-map-scattered-point?noredirect=1#comment245543_154975 — FelixIP, Oct 05 '15 at 18:20

Richard Fairhurst · Accepted Answer · 2015-10-05T18:31:51.033

Here is the code I have come up with, based originally on code posted by Curtis Price in post 13 of this thread: https://geonet.esri.com/thread/48226 The code reads a point's Floorplan field to find the corresponding floor plan polygon in a master feature class, then it duplicates the polygon into a new feature class, moves the polygon to the be centered on the point, and rotates the polygon based on a rotation angle field in the point.

You will need to modify the paths and feature class names in the Main section of the code and modify the sourceFieldList and fieldsList in the body of the code to fit your attribute names. This may not fit your needs, but it shows that it is possible to use points to control the location of polygons and should contain several methods that are suited to other applications that need to manipulate polygons.

# Name:         PolyToPtMoveRotate.py  
# Purpose:      Move centroid of polygon feature to point
#               and Rotate polygon based on the difference
#               between the polygon's rotation field
#               and the point's rotation field  
# Author:       Richard Fraihurst, rfairhur@rctlma.org  
# Created:      07/25/2015 09:44:25 AM  
# Environment:  ArcGIS 10.1+
# Credit:       Rotation based on RotateFC.py by Curtis Price
#               found in post 13 in this thread:
#               https://geonet.esri.com/thread/48226
# -------------------------------------------------------------------  

import os
import sys  
import traceback  
import arcpy
from arcpy import env  

def MoveRotatePolygonToPoint(inputPolygonLyr, inputPointLyr, outputFC):
    """Move and Rotate Polygon To Point  

    inputPolygonFC  Input polygon features  
    inputPointFC    Input point features  
    outputFC        Output feature class  

    As the output feature class no longer has a "real" xy locations,  
    after rotation, it no coordinate system defined.  
    """  

    def RotateXY(x, y, xc=0, yc=0, angle=0, units="DEGREES"):  
        """Rotate an xy cooordinate about a specified origin  

        x,y      xy coordinates  
        xc,yc   center of rotation  
        angle   angle  
        units    "DEGREES" (default) or "RADIANS"  
        """  
        import math  
        x = x - xc  
        y = y - yc  
        # make angle clockwise (like Rotate_management)  
        #angle = angle * -1  
        if units == "DEGREES":  
            angle = math.radians(angle)  
        xr = (x * math.cos(angle)) - (y * math.sin(angle)) + xc  
        yr = (x * math.sin(angle)) + (y * math.cos(angle)) + yc  
        return xr, yr  

    print("Starting method")
    # temp names for cleanup  
    env_file = None  
    lyrFC, lyrTmp, lyrOut   = [None] * 3  # layers  
    tmpFC  = None # temp dataset  
    Row, Rows, oRow, oRows = [None] * 4 # cursors  
    print("Temp variables set up")

    try:

        # verify input is correct shapetypes
        dPoly = arcpy.Describe(inputPolygonLyr)
        polyShpType = dPoly.shapeType
        if polyShpType != 'Polygon':
            raise Exception, "Shape type cannot be {0} for inputPolygonLyr".format(polyShpType) 

        dPoint = arcpy.Describe(inputPointLyr)
        ptShpType = dPoint.shapeType
        if ptShpType != 'Point':
            raise Exception, "Shape type cannot be {0} for inputPointLyr".format(ptShpType)
        FID = dPoint.OIDFieldName  

        if arcpy.Exists(outputFC):
            arcpy.Delete_management(outputFC)

        print("Deleted outputFC")

        # set up environment  
        env_file = arcpy.CreateScratchName("xxenv",".xml","file",  
                                           os.environ["TEMP"])  
        arcpy.SaveSettings(env_file)  

        # Disable any GP environment clips or project on the fly  
        #arcpy.ClearEnvironment("extent")  
        #arcpy.ClearEnvironment("outputCoordinateSystem")  

        WKS = env.workspace  
        if not WKS:  
            if os.path.dirname(outputFC):  
                WKS = os.path.dirname(outputFC)  
            else:  
                WKS = os.path.dirname(  
                    arcpy.Describe(inputPolygonFC).catalogPath)  
        env.workspace = env.scratchWorkspace = WKS  

        print("Set up scratch worksapce")

        # create temp feature class  
        tmpFC = arcpy.CreateScratchName("xxfc","","featureclass")  
        arcpy.CreateFeatureclass_management(os.path.dirname(tmpFC),  
                                            os.path.basename(tmpFC),  
                                            polyShpType)  
        lyrTmp = "lyrTmp"  
        arcpy.MakeFeatureLayer_management(tmpFC, lyrTmp)  
        # set up id field (used to join later)  
        TFID = "XXXX_FID"  
        arcpy.AddField_management(lyrTmp, TFID, "LONG")  
        arcpy.DeleteField_management(lyrTmp, "ID")  

        print("Created temp layer")

        sourceFieldsList = ['FLOORPLAN','SHAPE@','BUILDING_X_COORDINATE',
                                      'BUILDING_Y_COORDINATE','ROTATION']  

        # Use list comprehension to build a dictionary from a da SearchCursor  
        valueDict = {r[0]:(r[1:]) for r in arcpy.da.SearchCursor(inputPolygonLyr, sourceFieldsList)}  

        print("Read Polygon Template dictionary")

        fieldsList = ['SHAPE@','OID@','FLOORPLAN','BUILDING_X_COORDINATE',
                                      'BUILDING_Y_COORDINATE','ROTATION']
        # open read and write cursors  
        Rows = arcpy.da.SearchCursor(inputPointLyr,
                                     fieldsList)  

        print("Opened search cursor on points")

        insertFieldsList = ['SHAPE@',TFID]
        oRows = arcpy.da.InsertCursor(lyrTmp,insertFieldsList)  

        print("Opened Insert Cursors on output")

        parts = arcpy.Array()  
        rings = arcpy.Array()  
        ring = arcpy.Array()
        cnt = 0              
        print("Starting to create moved and rotated polygons")
        for Row in Rows:  
            #print(str(Row))
            floorplan = str(Row[2])
            if floorplan in valueDict:
                #print("Found {} in dictionary".format(floorplan))
                shp = valueDict[floorplan][0]
                #print(str(shp))
                xCen = float(Row[3])
                #print(str(xCen))
                xOffset = xCen - float(valueDict[floorplan][1])
                #print(str(xOffset))
                yCen = float(Row[4])
                #print(str(yCen))
                yOffset = yCen - float(valueDict[floorplan][2])
                #print(str(yOffset))
                angle = float(Row[5]) - float(valueDict[floorplan][3])
##                if angle < 0:
##                    angle = angle - 180
                #print(str(angle))
                #print("Read Fields")
                p = 0  
                for part in shp:  
                    for pnt in part:  
                        if pnt:
                            xMoved = pnt.X + xOffset
                            yMoved = pnt.Y + yOffset
                            x, y = RotateXY(xMoved, yMoved, xCen, yCen, angle)  
                            ring.add(arcpy.Point(x, y, pnt.ID))  
                        else:  
                            # if we have a ring, save it  
                            if len(ring) > 0:  
                                rings.add(ring)  
                                ring.removeAll()  
                    # we have our last ring, add it  
                    rings.add(ring)  
                    ring.removeAll()  
                    # if only one, remove nesting  
                    if len(rings) == 1: rings = rings.getObject(0)  
                    parts.add(rings)  
                    rings.removeAll()  
                    p += 1  

                # if only one, remove nesting  
                if len(parts) == 1: parts = parts.getObject(0)  
                if dPoly.shapeType == "Polyline":  
                    shp = arcpy.Polyline(parts)  
                else:  
                    shp = arcpy.Polygon(parts)  
                parts.removeAll()
                #print("Got new Shape")
                oRow = [shp, Row[1]]
                oRows.insertRow(oRow)
                cnt += 1

        del oRow, oRows # close write cursor (ensure buffer written)  
        oRow, oRows = None, None # restore variables for cleanup  

        print("Created {0} polygons".format(cnt))

        # join attributes, and copy to output  
        arcpy.AddJoin_management(lyrTmp, TFID, inputPointLyr, FID)
        print("Added Join")
        env.qualifiedFieldNames = False  
        arcpy.Merge_management(lyrTmp, outputFC)  
        print("Merged")
        lyrOut = "lyrOut"  
        arcpy.MakeFeatureLayer_management(outputFC, lyrOut)  
        # drop temp fields 2,3 (TFID, FID)  
        fnames = [f.name for f in arcpy.ListFields(lyrOut)]  
        dropList = ";".join(fnames[2:4])  
        arcpy.DeleteField_management(lyrOut, dropList)  
        print("Joined Attributes")

##    except MsgError, xmsg:  
##        arcpy.AddError(str(xmsg))  
    except arcpy.ExecuteError:  
        tbinfo = traceback.format_tb(sys.exc_info()[2])[0]  
        arcpy.AddError(tbinfo.strip())  
        arcpy.AddError(arcpy.GetMessages())  
        numMsg = arcpy.GetMessageCount()  
        for i in range(0, numMsg):  
            arcpy.AddReturnMessage(i)  
    except Exception, xmsg:  
        tbinfo = traceback.format_tb(sys.exc_info()[2])[0]  
        arcpy.AddError(tbinfo + str(xmsg))  
    finally:  
        # reset environment  
        if env_file: arcpy.LoadSettings(env_file)  
        # Clean up temp files  
        for f in [lyrFC, lyrTmp, lyrOut, tmpFC, env_file]:  
            try:  
                if f: arcpy.Delete_management(f)  
            except:  
                pass  
        # delete cursors  
        try:  
            for c in [Row, Rows, oRow, oRows]: del c  
        except:  
            pass  


if __name__ == '__main__':  
    inputPolygonFC = r"L:\rfairhur\Layers\Building Outlines\Building Outlines.gdb\Building_Outline_Templates"
    inputPolygonLyr = "inputPolygonLyr"
    arcpy.MakeFeatureLayer_management(inputPolygonFC, inputPolygonLyr)
    arcpy.SelectLayerByAttribute_management(inputPolygonLyr, "NEW_SELECTION", " GARAGE_SIDE > ' ' ")
    inputPointFC = r"L:\rfairhur\Layers\Building Outlines\Building Outlines.gdb\Building_Points"
    inputPointLyr = "inputPointLyr"
    arcpy.MakeFeatureLayer_management(inputPointFC, inputPointLyr)
    arcpy.SelectLayerByAttribute_management(inputPointLyr, "NEW_SELECTION", " GARAGE_SIDE > ' ' ")
    outputFC = r"L:\rfairhur\Layers\Building Outlines\Building Outlines.gdb\Building_Outline_Polygons_FC"
    MoveRotatePolygonToPoint(inputPolygonLyr, inputPointLyr, outputFC)

Is centering existing polygons on a point layer possible?

1 Answers1