4

I am looking for any guidance on how to calculate setbacks for a few parcels.

I am using ArcGIS 10.3

I have a shapefile of a few land parcel. They have columns of Side_Setback, Front_Setback, and Rear_Setback.

The image below explains setbacks well. enter image description here

I am looking to take the original parcel layer and have a model builder and/or python script cut those polygons into a polygon of that parcels buildable area and have the same attributes of the original polygon.

3
  • 1
    You probably need to add more of a description. Start with defining what a set back is. I assume you have legal land parcels? Maybe you could screen dump and annotate an example.
    – GISI
    Jan 20, 2016 at 4:45
  • Have you reviewed the Q&As returned by this search? gis.stackexchange.com/search?q=setback+parcel+is%3Aquestion
    – PolyGeo
    Jan 20, 2016 at 4:55
  • Those Q&As dont have any answers to them that would help solve my issue
    – GIS Pat
    Feb 9, 2016 at 1:20

1 Answer 1

2

Ok, if the parcels are all rectangles with one road frontage you could use a python script at a very low level to generate the setback lines and create a new polygon.

To do this you need to use geometry objects which give you the ability to manipulate parts of a feature and do various spatial operations between the parts.

I don't like recommending a loop around each parcel, but you say you haven't got many so it won't take too long. Here is my off the cuff pseudocode:

For each polygon:

  • collect attributes from polygon as a dictionary of setbacks
  • convert feature to a polygon geometry
  • split into a list of sides
  • identify each side (front, back, left, right) by using distanceTo of the boundary midpoints to road geometry objects
  • buffer each boundary by the offset
  • clip off the margin from the parcel geometry using each buffer
  • output a new polygon feature, adding back original key with
  • extra attributes

Here is a working script, 0.3 seconds per parcel on my laptop

#-------------------------------------------------------------------------------
# Name:        setback.py
# Purpose:     prototype setback calculation
#
# Author:      kimo
#
# Created:     09/02/2016
# Copyright:   (c) kimo 2016
# Licence:     Creative Commons 3.0 New Zealand
#-------------------------------------------------------------------------------
# only requires an ArcView licence

import sys
import arcpy
import datetime

try:
    parcel = sys.argv[1]
    road = sys.argv[2]
except:
    parcel_fc = 'd:/project/stackExchange/cadastral.gdb/parcel'
    road_fc = 'd:/project/stackExchange/cadastral.gdb/road'
debug = False
ws = 'd:/project/stackExchange/cadastral.gdb'
arcpy.env.workspace = ws
arcpy.env.overwriteOutput = True
start = datetime.datetime.now()
# get a list of road geometries to find nearest
g = arcpy.Geometry()
g_roads = arcpy.management.CopyFeatures(road_fc,g)
out_fc = 'build'
# open an insert cursor for output
if arcpy.Exists(out_fc):
    arcpy.management.Delete(out_fc)
arcpy.management.CreateFeatureclass(ws,'build','polygon',
              template=parcel_fc,spatial_reference=parcel_fc)
icur = arcpy.da.InsertCursor(out_fc,['SHAPE@','par_id','front_setback','rear_setback','side_setback'])
count = 0
with arcpy.da.SearchCursor(parcel_fc,['SHAPE@','par_id','front_setback','rear_setback','side_setback']) as cur:
    for row in cur:
        g_parcel = row[0]
        par_id = row[1]
        front = row[2]
        rear = row[3]
        side = row[4]
        if debug: print "parcel ID: {} front: {} rear: {} sides: {}".format(*row[1:])
        # get nearest road segment to polygon
        shortest = 1E6
        for r in g_roads:
            dist = r.distanceTo(g_parcel)
            if dist < shortest:
                r_seg = r
                shortest = dist
        if debug: print "Closest road segment is {:.2f} metres".format(shortest)

        bound = g_parcel.boundary()
        # only process exactly 5 points for rectangles
        if bound.pointCount == 5:
            pnt_array = bound.getPart(0) # only one part
            if debug: print 'points in polygon', len(pnt_array)
            # create a list of polyline geometry (with a spatial ref) sides
            # from an array of the end points for each side
            sides = []
            sr = arcpy.Describe(parcel_fc).spatialReference
            for n in range(4):
                sides.append(arcpy.Geometry('polyline',
                arcpy.Array([pnt_array[n],pnt_array[n+1]]),sr))
            # find midpoints for front (closest) and back (farthest), remainder are sides
            d = {}
            for s in sides:
                midpoint = arcpy.Geometry('point',s.centroid,sr)
                d[sides.index(s)] = midpoint.distanceTo(r_seg)
            if debug: print d
            order = sorted(d, key=d.get)
            if debug: print "road {} sides {}, {} back {}".format(*order)
            # now buffer each side by margin and subtract from original parcel poly
            front_buf = sides[order[0]].buffer(front)
            rear_buf = sides[order[3]].buffer(rear)
            side1_buf = sides[order[1]].buffer(side)
            side2_buf = sides[order[2]].buffer(side)
            p1 = g_parcel.difference(front_buf)
            p2 = p1.difference(rear_buf)
            p3 = p2.difference(side1_buf)
            p4 = p3.difference(side2_buf).projectAs(sr)
            # export using an insert cursor with attributes
            outrow = [p4,par_id,front,rear,side]
            icur.insertRow(outrow)
            count+=1
del cur
del icur
print "Finished {} in  {}".format(count, datetime.datetime.now() - start)

output showing id and offsets



#-------------------------------------------------------------------------------
# Name:        setback_n.py
# Purpose:     prototype setback calculation
#              generalised for more than 4 sides
# Author:      kimo
#
# Created:     15/02/2016
# Copyright:   (c) kimo 2016
# Licence:     Creative Commons 3.0 New Zealand
#-------------------------------------------------------------------------------
# only requires an ArcView licence

import sys
import arcpy
import datetime

try:
    parcel = sys.argv[1]
    road = sys.argv[2]
except:
    parcel_fc = 'd:/project/stackExchange/cadastral.gdb/parcel'
    road_fc = 'd:/project/stackExchange/cadastral.gdb/road'
debug = False
ws = 'd:/project/stackExchange/cadastral.gdb'
arcpy.env.workspace = ws
arcpy.env.overwriteOutput = True
start = datetime.datetime.now()
# get a list of road geometries to find nearest
g = arcpy.Geometry()
g_roads = arcpy.management.CopyFeatures(road_fc,g)
out_fc = 'build'
# open an insert cursor for output
if arcpy.Exists(out_fc):
    arcpy.management.Delete(out_fc)
arcpy.management.CreateFeatureclass(ws,'build','polygon',template=parcel_fc,spatial_reference=parcel_fc)
icur = arcpy.da.InsertCursor(out_fc,['SHAPE@','par_id','front_setback','rear_setback','side_setback'])
count = 0
with arcpy.da.SearchCursor(parcel_fc,['SHAPE@','par_id','front_setback','rear_setback','side_setback']) as cur:
    for row in cur:
        g_parcel = row[0]
        par_id = row[1]
        front = row[2]
        rear = row[3]
        side = row[4]
        if debug: print "parcel ID: {} front: {} rear: {} sides: {}".format(*row[1:])
        # get nearest road segment to polygon
        shortest = 1E6
        for r in g_roads:
            dist = r.distanceTo(g_parcel)
            if dist = 5:
            pnt_array = bound.getPart(0) # only one part
            if debug: print 'points in polygon', len(pnt_array)
            # create a list of polyline geometry (with a spatial ref) sides
            # from an array of the end points for each side
            sides = []
            sr = arcpy.Describe(parcel_fc).spatialReference
            vertices = bound.pointCount
            if debug: print "Vertices", vertices
            for n in range(vertices-1):
                sides.append(arcpy.Geometry('polyline',
                arcpy.Array([pnt_array[n],pnt_array[n+1]]),sr))
            # find midpoints for front (closest) and back (farthest), remainder are sides
            d = {}
            side_buf = {}
            for s in sides:
                midpoint = arcpy.Geometry('point',s.centroid,sr)
                d[sides.index(s)] = midpoint.distanceTo(r_seg)
            if debug: print d
            order = sorted(d, key=d.get)
            if debug: print "road {} sides {}, {} back {}".format(*order)
            # now buffer each side by margin and subtract from original parcel poly
            front_buf = sides[order[0]].buffer(front)
            rear_buf = sides[order[-1]].buffer(rear)
            for n in range(1,vertices-2):
                side_buf[n] = sides[order[n]].buffer(side)
            p = g_parcel.difference(front_buf)
            for n in range(1,vertices-2):
                p = p.difference(side_buf[n])
            p = p.difference(rear_buf).projectAs(sr)

            # export using an insert cursor with attributes
            outrow = [p,par_id,front,rear,side]
            icur.insertRow(outrow)
            count+=1
del cur
del icur
print "Finished {} in  {}".format(count, datetime.datetime.now() - start)


enhanced for more than four sides

Some irregular shaped parcels show the limits of this method to identify which boundaries are front and rear. Also the short neck does not cleanly draw a line extending across the right-of-way. Getting a bit tricky to handle all these edge cases with a set of coherent rules. It just goes to show how planners make impractical zoning rules!

If each parcel has different rules, depending on the zoning and landuse then I suppose you have to do it per parcel. If you have complex boundaries you will need to identify the boundary line types. This used to be easy with coverages because there was a PAL polygon-arc-list that kept a relation between polygons and boundaries. The only way with arcgis is to use two layers with some topology rules to manage the duplication of polylines and polygons to warn when things are not in sync. If side boundaries have different offsets you will need two coincident polylines and a tag to identify the correct parcel. Yuck.

2
  • Would there possibly be any work around this as the vast majority of my parcels are not rectangles
    – GIS Pat
    Feb 11, 2016 at 3:24
  • Provided that there is a clear front boundary and rear boundary, all the other boundaries would be treated as side boundaries. None of the angles need to be right angles because the line creates a buffer. Back properties with a narrow entrance would be ok. I would need some examples that were not Ok to see what could be done. The main requirement is that there are not irregular boundaries with lots of vertices, they have to be Right Lines (in survey jargon)
    – kimo
    Feb 15, 2016 at 8:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.