python select features just interiors of geometry intersect

Question

Hi folks,

Does python provide funtion to select features jjust interiors of geometry intersect? Similar like ArcOject esriSeptialRelationEnum = 2 which Intersection excluding boundary touch.

I try to several ways and can't figure it out.

any idea?

thanks
choilon

KimOllivier · Answer

Yes, that is easy to select. You just have to do your own logic with the selected sets to add or subtract until you are left with the cases you are interested in. If you can logically describe to me what you need, you can create a program to do the same. That is what makes programming so fascinating.

First you select intersected features (INTERSECT), and from that subset (using a layer) you select features that are touching (BOUNDARY_TOUCHES).
You can then also remove contained polygons from the selection using CONTAINS.

In my case I made up a separate selection to identify a group that I did not want (touching but outside) and used an attribute selection to remove them from the intersecting and inside set.

It does not matter that an operator selects other cases you are not interested in if you have already selected them out so that you know that those cases do not apply for your set.
Try it interactively using the SelectbyLocation option on some samples to get the logic.

longcui · Answer

Hi Kimo,

I think the diagrams from A-M do not include the scenario like the two features intersect and share boundary, but not contain each other (like the image attached). In that case, when exclude DJI, this scenario is excluded as well.
[ATTACH=CONFIG]22359[/ATTACH]
thanks
Long

example.JPG

KimOllivier · Answer

I used the diagrams in the help called Select By Location: graphic examples.
These have nice summaries of each selection operator summaries in a set of diagrams A - M

Suppose you want to exclude polygons touching on the outside (diagram I)

Case I is covered in INTERSECT, BOUNDARY_TOUCHES, SHARE_A_LINE_SEGMENT_WITH
Other operators do not include diagram I so you add and remove using other operators until you are left with case I.
I made up a matrix of the letters so that I could see which ones to add or subtract.

Since we do not want to confuse the list with features that are not touching at all, I start with an intersect subset and make a temporary featureclass. At the end I simply get the remaining keys (par_id) of the features required and make a normal SQL selection using a "KEY IN (key,key,key)" expression.

This snippet of code illustrates the flaming hoops that you have to jump through compared to ARC/INFO coverage selections.
Since the new data is not topologically correct I have to add a snap as well to avoid errors from sliver overlaps.
# ...snip 
# get all intersections ACDEFGHIJKM (see help diagrams for letters)
arcpy.management.SelectLayerByLocation(fcLayer,'INTERSECT',gaLayer,selection_type='NEW_SELECTION')
# copy to a temp layer, snap, then remove those outside but touching
snapFC = prefix+fcAbbrev+"_snap"
arcpy.management.CopyFeatures(fcLayer,snapFC)
arcpy.management.MakeFeatureLayer(snapFC,touchLayer) # Have no selection to start
arcpy.management.SelectLayerByLocation(gaLayer,'INTERSECT',fcLayer,selection_type='ADD_TO_SELECTION')
arcpy.edit.Snap(touchLayer,[[gaLayer,'VERTEX',4],[gaLayer,'EDGE',4]])
# get outside touching to exclude DIJ (help diagram cases)
arcpy.management.SelectLayerByLocation(touchLayer,'BOUNDARY_TOUCHES',gaLayer,selection_type='NEW_SELECTION')
arcpy.management.SelectLayerByLocation(touchLayer,'WITHIN',gaLayer,selection_type='REMOVE_FROM_SELECTION')
arcpy.management.SelectLayerByLocation(touchLayer,'CONTAINS',gaLayer,selection_type='REMOVE_FROM_SELECTION')
# subtract outside touching if any DIJ (see help)
cTouch = int(arcpy.management.GetCount(touchLayer).getOutput(0))
if cTouch == 0:
    arcpy.analysis.SpatialJoin(fcLayer,gaLayer, prefix+fcAbbrev,'JOIN_ONE_TO_MANY','KEEP_COMMON',"#",'INTERSECT')
    count = int(arcpy.management.GetCount(prefix+fcAbbrev).getOutput(0))
    if count > 0:
        dCount[fcAbbrev] = count
        msg = "GA %s %s %s %d" % (aggregate_area,arcpy.management.GetCount(gaLayer),fcAbbrev,dCount[fcAbbrev])
        logging.info(msg)
        print msg
    else:
        arcpy.management.Delete(prefix+fcAbbrev)
else:
    lstPar = [row.par_id for row in arcpy.SearchCursor(touchLayer)]
    if len(lstPar) > 0:
        if len(lstPar) == 1:
            expr = "PAR_ID = "+str(lstPar[0])
        else:
            expr = "PAR_ID in "+str(tuple(lstPar))
        # print "touching outside",cTouch #,expr
        # remove from list of new parcels of interest
        try:
            arcpy.management.SelectLayerByAttribute(fcLayer,"REMOVE_FROM_SELECTION",expr)
        except:
            print "error",expr
    #
    if int(arcpy.management.GetCount(fcLayer).getOutput(0)) > 0:
        arcpy.analysis.SpatialJoin(fcLayer,gaLayer, prefix+fcAbbrev,'JOIN_ONE_TO_MANY','KEEP_COMMON',"#",'INTERSECT')
        count = int(arcpy.management.GetCount(prefix+fcAbbrev).getOutput(0))
        if count > 0:
            dCount[fcAbbrev] = count
            # tag touched GA layers
            arcpy.management.SelectLayerByLocation(gaLayer,'INTERSECT',fcLayer,selection_type='ADD_TO_SELECTION')
            msg = "GA %s %s %s %d" % (aggregate_area,arcpy.management.GetCount(gaLayer),fcAbbrev,dCount[fcAbbrev])
            logging.info(msg)
            print msg
        else:
            arcpy.management.Delete(prefix+fcAbbrev)