1

I've tried digging around for an answer to this question, but I've only been able to find answers that are either too broad or too different from my case, so here goes!

Context:

I have a layer of city parcels, and I am trying to find all the smaller parcels that lie within bigger container parcels (i.e. I am trying to identify condos, trailers, or any other type of housing that has its own unique address, owner name, and ID that is different from the container parcel). The image below illustrates my task: the red parcel is the container and the blue circled parcel is one of several that I am trying to programmatically identify.

Completely enclosed parcels

My goal is to create a dictionary where the key of each entry is the smaller parcel's ID# and the value is the container parcel's ID# (basically a relational table of bounded parcels).

I have tried working with the Polygon Neighbors tool, but it's redundant in my case for two reasons:

  1. At the city-scale, there are ~44,000 parcels, and each of those parcels has anywhere from one to tens of neighbors, and the list would be frighteningly endless. You would think that I could just run the tool and then query out all the entries that only have one neighbor (meaning they only have one enclosing neighbor), but...
  2. There are small parcels within larger parcels that are not completely enclosed AND have more than one neighbor (see picture below), so even if I ran the tool and weeded out table keys with repeating IDs, I would risk losing some info.

Not-completely enclosed parcels

How I think the problem can be solved:

My fundamental idea has been to compare extents between and among parcels. Essentially, I'm thinking that if the bounding box of the smaller parcel is completely within the bounding box of the bigger parcel, I can have my script write those IDs to a dictionary. As far as coding that workflow, it seems that I should make use of two Search Cursors on the same feature class and compare geometries line-by-line.

My code:

parcels = path_to_testing_fc

parcelDict = {}
with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs1, sCurs2:
    for inside in sCurs1: #small parcels' loop
        for container in sCurs2: #container parcels' loop
            if inside[0].extent.within(container[0], 'PROPER'): #"proper" ensures my completely within rule
                parcelDict[str(inside[1])] = 'Contains ' + str(container[1])

When I run this, I get no errors, but the "parcelDict" dictionary is also empty, so I'm kind of in limbo.

Is the double search cursor a viable, efficient way to get this task done?In the same vein, is there a better, more efficient way to accomplish this goal?

Getting the snippet to work on my small test feature class is entirely different from getting the code to run through 44,000 parcels in a timely manner. Perhaps I'll look at the code for Polygon Neighbors and see if I can't glean something.


I opened a new console and my code from above suddenly started throwing errors (I probably forgot to delete cursors from other attempts at solving this problem), so here's an update on what my code looks like (still not giving errors and still not populating my dictionary...)

parcels = path_to_test_fc

parcelDict = {}
with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs1:
    for inside in sCurs1: #small parcels' loop
        with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs2:
            for container in sCurs2: #container parcels' loop
                if inside[0].extent.within(container[0], 'PROPER'):
                    parcelDict[str(inside[1])] = 'Contains ' + str(container[1])
2

Here is the code snippet that finally solved my problem:

with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs1:
    for inside in sCurs1: #small parcels' loop
        xmax1 = inside[0].extent.XMax; ymax1 = inside[0].extent.YMax
        xmin1 = inside[0].extent.XMin; ymin1 = inside[0].extent.YMin
        with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs2:
            for container in sCurs2: #container parcels' loop
                xmax2 = container[0].extent.XMax; ymax2 = container[0].extent.YMax
                xmin2 = container[0].extent.XMin; ymin2 = container[0].extent.YMin
                    if xmax1 < xmax2 and xmin1 > xmin2 and ymax1 < ymax2 and ymin1 > ymin2:
                        parcelDict[str(inside[1])] = 'Is inside ' + str(container[1])

The reason this works as opposed to my original code is because the container parcels I am interested in have holes where the smaller parcels are located. The new code now compares geometries rather than using ESRI's extent.within, since those holes throw that test off.


Here is the code snippet that proved most efficient:

parcels = path_to_parcel_in_sde
parcelDict = {}

with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs:
    for row in sCurs:
        xmax = row[0].extent.XMax; ymax = row[0].extent.YMax
        xmin = row[0].extent.XMin; ymin = row[0].extent.YMin
        parcelDict[str(row[1])] = [xmax, ymax, xmin, ymin]

finalDict = {}
for key1, val1 in parcelDict.iteritems():
    setDict = {k:v for k,v in parcelDict.iteritems() if k[:-3] == key1[:-3]}
    for key2, val2 in setDict.iteritems():
        if val1[0] < val2[0] and val1[1] < val2[1] and val1[2] > val2[2] and val1[3] > val2[3]:
            finalDict[key1] = 'Is inside ' + key2

It uses one Search Cursor to create a dictionary, and then uses dictionary comprehension to create a subset of parcels that have similar PARCEL_ID values to compare. On my test feature class, the first snippet took 18.4 seconds, while this snippet took 0.4 seconds. I used this snippet on the entire parcel layer and it took 6 minutes to process all 44,000 parcels.

0

I am assuming that all parcels are within the same feature class.

To solve this I would:

  1. Use the Polygon Neighbors tool to create a table of the neighboring parcels
  2. Use dictionary comprehension to build a dictionary of each parcel with a list of its neighbors
  3. Use that dictionary as the basis for completing any other tasks that need to know a polygon's neighbors.
  • it's a start, but the problem I run into is that this solution gives me a list of all parcels and their neighbors, when I just need ones that contain other parcels. By my estimates, there are about 1000 parcels with others inside them, so a dictionary of 1000 would be preferable to one of 44000. – Jesse Nestler Feb 7 '18 at 14:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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