Given multiple polygons that overlap in multiple ways, I would like to export from these features all polygons that don't overlap with others, iteratively.
The product would be a number of features with no overlap that, when summed together, make up the original.
The products could then be used as input to Zonal Statistics, and this would be much faster than iterating Zonal Statistics over each polygon.
I have been trying to code this in ArcPy without success.
Does code to do this already exist?
Answer
The methodology recommended by #whuber inspired me to take a new direction, and here is my arcpy solution, in two functions. The first, called countOverlaps, make two fields, "overlaps" and "ovlpCount" to record for each poly which polys overlapped with it, and how many overlaps occurred. The second function, explodeOverlaps, creates a third field, "expl", which gives a unique integer to each group of non-overlapping polys. The user can then export new fc's based on this field. The process is broken into two functions because I think the countOverlaps tool can prove useful by itself. Please excuse the sloppiness of the code (and the careless naming convention), as it's pretty preliminary, but it works. Also make sure that the "idName" field represents a field with unique IDs (only tested with integer IDs). Thank you whuber for providing me with the framework necessary to approach this problem!
def countOverlaps(fc,idName):
intersect = arcpy.Intersect_analysis(fc,'intersect')
findID = arcpy.FindIdentical_management(intersect,"explFindID","Shape")
arcpy.MakeFeatureLayer_management(intersect,"intlyr")
arcpy.AddJoin_management("intlyr",arcpy.Describe("intlyr").OIDfieldName,findID,"IN_FID","KEEP_ALL")
segIDs = {}
featseqName = "explFindID.FEAT_SEQ"
idNewName = "intersect."+idName
for row in arcpy.SearchCursor("intlyr"):
idVal = row.getValue(idNewName)
featseqVal = row.getValue(featseqName)
segIDs[featseqVal] = []
for row in arcpy.SearchCursor("intlyr"):
idVal = row.getValue(idNewName)
featseqVal = row.getValue(featseqName)
segIDs[featseqVal].append(idVal)
segIDs2 = {}
for row in arcpy.SearchCursor("intlyr"):
idVal = row.getValue(idNewName)
segIDs2[idVal] = []
for x,y in segIDs.iteritems():
for segID in y:
segIDs2[segID].extend([k for k in y if k != segID])
for x,y in segIDs2.iteritems():
segIDs2[x] = list(set(y))
arcpy.RemoveJoin_management("intlyr",arcpy.Describe(findID).name)
if 'overlaps' not in [k.name for k in arcpy.ListFields(fc)]:
arcpy.AddField_management(fc,'overlaps',"TEXT")
if 'ovlpCount' not in [k.name for k in arcpy.ListFields(fc)]:
arcpy.AddField_management(fc,'ovlpCount',"SHORT")
urows = arcpy.UpdateCursor(fc)
for urow in urows:
idVal = urow.getValue(idName)
if segIDs2.get(idVal):
urow.overlaps = str(segIDs2[idVal]).strip('[]')
urow.ovlpCount = len(segIDs2[idVal])
urows.updateRow(urow)
def explodeOverlaps(fc,idName):
countOverlaps(fc,idName)
arcpy.AddField_management(fc,'expl',"SHORT")
urows = arcpy.UpdateCursor(fc,'"overlaps" IS NULL')
for urow in urows:
urow.expl = 1
urows.updateRow(urow)
i=1
lyr = arcpy.MakeFeatureLayer_management(fc)
while int(arcpy.GetCount_management(arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",'"expl" IS NULL')).getOutput(0)) > 0:
ovList=[]
urows = arcpy.UpdateCursor(fc,'"expl" IS NULL','','','ovlpCount D')
for urow in urows:
ovVal = urow.overlaps
idVal = urow.getValue(idName)
intList = ovVal.replace(' ','').split(',')
for x in intList:
intList[intList.index(x)] = int(x)
if idVal not in ovList:
urow.expl = i
urows.updateRow(urow)
ovList.extend(intList)
i+=1
No comments:
Post a Comment