Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this community
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.
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.
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:
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)
#-------------------------------------------------------------------------------
# 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)
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.
