How to identify donut holes with the Field Calculator?

[AndrewRudin]

Anyone have any ideas why a field calculator function would return different results than an arcpy function written in pyscripter?

I'm processing hundreds of thousands of polygons in ArcGIS Desktop 10.2.2 in a model and I want to be able to capture which ones contain a donut hole (a.k.a. at least one part of the polygon forms at least one ring).  Since I don't see any pre-built tools for this, I assume the only way to do this without ArcObjects is through an iteration of each record's shape geometry using Python.  I built a python function in pyscripter using an example in the ESRI help at this location as a template and it worked on my sample data.

My problem now is that I get different results when running the python function in ArcMap's field calculator. When the code reads a feature that does have a donut hole in it, the code says it doesn't.

Here's my PyScripter code that gives correct results:

import arcpy

#Begin donut finder function

def has_donut(feat):
    donut = 'F'
    partcount = feat.partCount
    for p in range (0,partcount):
        part = feat.getPart(p)
        for pnt in part:
            if not pnt:
                donut = 'T'
                return donut
    return donut

#set input polygon feature class

in_fc = r"C:\my_geodatababase.gdb\my_feature_class

#Loop over feature class records and print donut search result

with arcpy.da.SearchCursor(in_fc,['OID@','SHAPE@']) as cur:
    for row in cur:
        oid = str(row[0])
        donut_type = has_donut(row[1])
        print oid + ": " + donut_type

Here's my ArcMap field calculator setup that returns False even when a polygon contains a donut:

Parser: Python

Pre-Logic Script Code:

def has_donut(feat):

    donut = 'F'

    partcount = feat.partCount

    for p in range (0,partcount):

        part = feat.getPart(p)

        for pnt in part:

            if not pnt:

                donut = 'T'

                return donut

    return donut

Field Calculator Expression:

has_donut( !SHAPE!)

[BruceHarold]

Hi

The Data Interoperability extension Workbench application has a HoleCounter transformer for this.

It will return 0 if a polygon has no holes, or the count of them if it does.

Regards

[DarrenWiens2]

I'm surprised this works in your original script:

for pnt in part:
    if not pnt:

For every time pnt exists in part, if it doesn't exist, do something.

[XanderBakker]

I can reproduce the problem, but I don't understand why the function behaves differently in the FieldCalculator.

I would suggest to use the script you build in PyScripter include parameters and create a tool from it and include it in your model... or better yet convert the model to python, include your logic and run from PyScripter.

I hardly ever use the field calculator (must be for something).

[AndrewRudin]

Thanks for the recommendation, Xander.  I think a toolbox exposed python script is probably the best way to go for my situation.  I wish the field calculator had a debugging function that could call pyscripter for issues like this.  My only guess on the difference is that the way field calculator loops over records is different than the searchcursor in arcpy, and thus the variables are getting handled differently.

Darren, on the issue of the "if not pnt" line and how it works to capture donut holes, I ran the loop again and set some watches:

Here's the watch list while its on a valid point in the part:

And here's what it looks like when it gets to the donut hole:

The only reason I knew about it is from the ESRI Help example about reading polygon geometries.  They had that if statement in there with a comment that it specifically is for inner rings.  It still begs the question what the inner ring really is, but I'm not that concerned with it at this juncture in my life.

Thank you everyone for your help

[DanPatterson_Retired]

A bit of a rant...I don't need explanations...I have heard them all.   However, I still marvel by the consistency with which things are instantiated.  For example, there is only one in the list that makes sense to me.

>>> # when is nothing really nothing
>>>
>>> pnt = arcpy.Point()   # a null point...or one would think
>>> arr = arcpy.Array()   # ditto for an array object
>>> xy = arcpy.Point(10,10)  # the real point
>>> nun = None
>>> lst = []
>>> # here we go...
>>> arr.extend([pnt,pnt,xy,nun,lst,[pnt,pnt,pnt]])
>>> arr
<Array [<Point (0.0, 0.0, #, #)>, <Point (0.0, 0.0, #, #)>, <Point (10.0, 10.0, #, #)>, None, <Array []>, <Array [<Point (0.0, 0.0, #, #)>, <Point (0.0, 0.0, #, #)>, <Point (0.0, 0.0, #, #)>]>]>
>>>
>>> for i in range(arr.count):
...  val = arr.getObject(i)
...  print('type and value: {}: {}'.format(type(val),val))
... 
type and value: <class 'arcpy.arcobjects.arcobjects.Point'>: 0 0 NaN NaN
type and value: <class 'arcpy.arcobjects.arcobjects.Point'>: 0 0 NaN NaN
type and value: <class 'arcpy.arcobjects.arcobjects.Point'>: 10 10 NaN NaN
type and value: <type 'NoneType'>: None
type and value: <class 'arcpy.arcobjects.arcobjects.Array'>: <geoprocessing array object object at 0x0D80EC68>
type and value: <class 'arcpy.arcobjects.arcobjects.Array'>: <geoprocessing array object object at 0x0D80EA88>
>>>

I miss null points...ie a point object that has no coordinates... or even NaN would be better than what was presented in the previous.
Can't null points (pnt above) and null arrays be instantiated as such.  Even the null list gets cast to an array, which should be a null array and the null array full of null points is a living breathing entity full of nothingness.

Can this cause problems?

>>> minimum = min([arr.getObject(i).X for i in range(3)])
>>> minimum
0.0
>>>

Yes...the minimum X coordinate of 2 null points and one real point is the non-existent X coordinate of the null point.
Later...going back to work with numpy arrays....