Monday, 13 August 2018

arcgis desktop - Distinguishing lines that intersect from lines that touch?


How do I distinguish between these cases in ArcGIS 10?



  • Case 1: Both endpoints of a line touch another line

  • Case 2: Both endpoints dangle over the lines it intersects


I am looking at the Trim Line function but that's not what I want (destructive).


The real world use case is to distinguish between service roads connecting both roadways of a freeway, and other cases of roads intersecting with freeways.


enter image description here enter image description here



Answer




For a single feature at a time, you can do this pretty easily interactively using the normal Select By Location dialog, using the following key as a guide to the spatial relationship types for line on line overlays (from Select by Location: graphic examples):



image
(source: arcgis.com)


Select line using line



INTERSECT A, C, D, E, F, G, H, I, J

CONTAINS G, H


COMPLETELY_CONTAINS G

CONTAINS_CLEMENTINI G, H

WITHIN F, H

COMPLETELY_WITHIN F

WITHIN_CLEMENTINI F, H


ARE_IDENTICAL_TO H

BOUNDARY_TOUCHES C, E

The relevant relationship types in this case are INTERSECT and BOUNDARY_TOUCHES. As you can see from the diagram above, you can use BOUNDARY_TOUCHES to select the features that touch an endpoint of the line. If exactly two features are selected then you have your Case 1. If a feature is not touched by any other features but only intersected by them, then BOUNDARY_TOUCHES will select nothing. INTERSECT will select all features that intersect regardless of whether they touch at an endpoint or not. So if you know there are no features touching endpoints, but you find there are features intersecting, then you have your Case 2.


To automate the process you can use the following Python script (implement as a script tool if desired) to calculate the number of touches and intersections for each feature in a feature class or layer:


import arcpy

################################ Configuration #################################
numTouchesField = "NUM_TOUCHES"
numIntersectionsField = "NUM_INTERSECTIONS"

################################################################################

def countTouches(layer, feature):
"""Returns the number of times the boundary of a feature touches other
features in the same feature layer."""
return countSpatialRelation(layer, feature, "BOUNDARY_TOUCHES")

def countIntersections(layer, feature):
"""Returns the number of times a feature intersects other features in the
same feature layer."""

return countSpatialRelation(layer, feature, "INTERSECT") - 1 # Subtract 1 because the feature will always intersect its clone in the feature layer

def countSpatialRelation(layer, feature, relation):
"""Returns the number of times a feature meets the specified spatial
relationship with other features in the same feature layer."""
arcpy.SelectLayerByLocation_management(layer, relation, feature)
count = int(arcpy.GetCount_management(layer).getOutput(0))
return count

def addField(table, fieldName, fieldType):

"""Adds a fields of the given name and type to a table, unless a field with
the same name already exists."""
desc = arcpy.Describe(table)
fieldInfo = desc.fieldInfo
fieldIndex = fieldInfo.findFieldByName(fieldName)
if fieldIndex == -1:
# Field does not exist, add it
arcpy.AddField_management(table, fieldName, fieldType)

def countTouchesAndIntersections(layer):

"""Adds and populates fields describing the number of times each feature
touches and intersects other features in the feature layer."""
addField(layer, numTouchesField, "LONG")
addField(layer, numIntersectionsField, "LONG")
desc = arcpy.Describe(layer)
shapeField = desc.shapeFieldName
rows = arcpy.UpdateCursor(layer)
for row in rows:
feature = row.getValue(shapeField)
row.setValue(numTouchesField, countTouches(layer, feature))

row.setValue(numIntersectionsField, countIntersections(layer, feature))
rows.updateRow(row)
del row, rows

if __name__ == "__main__":
layer = arcpy.MakeFeatureLayer_management(arcpy.GetParameterAsText(0))
countTouchesAndIntersections(layer)

Once that has run, you can easily query for the features that touch exactly twice and intersect exactly twice (Case 1), and those that touch 0 times and intersect exactly twice (Case 2).


Example definition queries:




  • Case 1 (Touches twice, intersects twice): "NUM_TOUCHES" = 2 AND "NUM_INTERSECTIONS" = 2

  • Case 2 (Touches none, intersects twice): "NUM_TOUCHES" = 0 AND "NUM_INTERSECTIONS" = 2


See the below screenshot for an illustration of the instances of the two cases being found: ArcMap screenshot showing various line intersection/touch relationships


Note that with real world data, normally street segments are broken up at intersections, and dangles only occur when roads pass over one another as at an interchange or bridge. So normally you have the same number of features intersecting as touching.


For the more general case, you might want to look for any dangles by checking whether "NUM_INTERSECTIONS" > "NUM_TOUCHES".


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...