Spacial Join closest line with a certain angle. Is it possible?

2330
6
12-20-2016 08:38 AM
JúlioLuta
New Contributor III

Hi,

Is it possible to spacial join a line feature set to another line feature set by proximity (closest) AND angle (ignoring all lines that have an angle/direction difference higher than X)?

From what I read about the Spacial Join Tool it seems it can only do the first, but wanted to check with those with more ArcGIS experience.

If it's indeed not possible, I 'm having doubts what the best solution would be:

1- Join each line with all that are within a certain distance (keeping the angle field of all) and then calculate which one has the most similar angle - but since I want to know the closest line that's in a certain angle interval in relation to the target line, I think this could give me the wrong results (especially as the distance radius increases)

2- Make and merge several Spacial Joins (closest) filtering the target feature ( angle=n ) and join feature ( angle=n+-x ) and just chose a x and go from n=0 to n=359? But I don't how to code this.

Cheers

0 Kudos
6 Replies
ClintonDow1
Occasional Contributor II

Hello Júlio,

You are correct that the spatial join tool doesn't have this functionality built in. I think an approach like you have in your first possible solution would work. Get the angle of a line using the function below:

import numpy as np


def angle_between(p1, p2):
    ang1 = np.arctan2(*p1[::-1])
    ang2 = np.arctan2(*p2[::-1])
    return np.rad2deg((ang1 - ang2) % (2 * np.pi))

Create a list of other lines which satisfy the spatial join. Then do a loop through the second list and use the above function to compare the angles of the lines. 

DanPatterson_Retired
MVP Emeritus

Another incarnation which scales well in terms of using real projected coordinates rather than those centered about an origin, is given below.  The difference is that you can specify an 'origin' and a 'destination' and the differences in the coordinates are used rather than the actual coordinates themselves.

You can use the following to either determine the angle that a line is directed towards relative to the x-axis or as an azimuth relative to north.  In this case, the origin and destination would be the start and end of the line.  

Alternately, you can pick a point from an origin feature and get the direction to one of a destination features points.  The attached would have be incorporated into a larger script to select the origin point from the origin features and the same for the destination feature.  You get the idea I hope (kind of like the Near tool, which could also be swung into service if you have the appropriate license level.

Anyway, here is the code followed by some results for a 0,0 origin and points located at compass points within root(2) of the origin.  The 'demo' function is written verbosely so you can see the construction.

def line_dir(orig, dest, fromNorth=False):
    """Direction of a line given 2 points
    : orig, dest - two points representing the start and end of a line.
    : fromNorth - True or False gives angle relative to x-axis)
    :Notes:
    :
    """
    orig = np.asarray(orig)
    dest = np.asarray(dest)
    dx, dy = dest - orig
    ang = np.degrees(np.arctan2(dy, dx))
    if fromNorth:
        ang = np.mod((450.0 - ang), 360.)
    return ang

def demo(xc=0, yc=0, fromNorth=True):
    """ run the demo with the data below """
    p0 = np.array([xc, yc])  # origin point
    p1 = p0 + [-1, 1]  # NW
    p2 = p0 + [0, 1]   # N
    p3 = p0 + [1, 1]   # NE
    p4 = p0 + [1, 0]   # E
    p5 = p0 + [1, -1]  # SE
    p6 = p0 + [0, -1]  # S
    p7 = p0 + [-1, -1] # SW
    p8 = p0 + [-1, 0]  # W
    #
    od = [[p0, p1], [p0, p2], [p0, p3], [p0, p4],
          [p0, p5], [p0, p6], [p0, p7], [p0, p8]]
    for pair in od:
        orig, dest = pair
        ang = line_dir(orig, dest, fromNorth=fromNorth)
        if fromNorth:
            dir = "From N."
        else: 
            dir = "From x-axis"
        args = [orig, dest, dir, ang]
        print("orig: {}: dest: {!s:<8} {}: {!s:>6}".format(*args))

# ---------------------------------------------------------------------
if __name__ == "__main__":
    """Main section...   """
    # print("Script... {}".format(script))
    xc = 300000   # pick an origin x  try:   0 or 300000 for example
    yc = 5025000  # pick an origin y  ditto: 0 or 5025000
    demo(xc, yc, fromNorth = True)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Some results with using the two sample origins.

Note! This can only be used for projected coordinates.

orig: [0 0]: dest: [-1  1]  From N.:  315.0
orig: [0 0]: dest: [0 1]    From N.:    0.0
orig: [0 0]: dest: [1 1]    From N.:   45.0
orig: [0 0]: dest: [1 0]    From N.:   90.0
orig: [0 0]: dest: [ 1 -1]  From N.:  135.0
orig: [0 0]: dest: [ 0 -1]  From N.:  180.0
orig: [0 0]: dest: [-1 -1]  From N.:  225.0
orig: [0 0]: dest: [-1  0]  From N.:  270.0

orig: [ 300000 5025000]: dest: [ 299999 5025001] From N.:  315.0
orig: [ 300000 5025000]: dest: [ 300000 5025001] From N.:    0.0
orig: [ 300000 5025000]: dest: [ 300001 5025001] From N.:   45.0
orig: [ 300000 5025000]: dest: [ 300001 5025000] From N.:   90.0
orig: [ 300000 5025000]: dest: [ 300001 5024999] From N.:  135.0
orig: [ 300000 5025000]: dest: [ 300000 5024999] From N.:  180.0
orig: [ 300000 5025000]: dest: [ 299999 5024999] From N.:  225.0
orig: [ 300000 5025000]: dest: [ 299999 5025000] From N.:  270.0
JúlioLuta
New Contributor III

Thanks Clinton and Dan!

As soon as I try the functions I'll come back with more detailed feedback.

Happy Holidays! 

0 Kudos
JúlioLuta
New Contributor III

Hi,

I wanted to share my progress in case it's useful to someone in the future. I also have a question because I can't understand why one step is not working.

I calculated the angle of the two feature classes using the following code I found online and which uses the !shape! and the field calculator:

# Pre-Loci Script Code
import math
def GetGeographicalDegrees(shape):
  radian = math.atan2(shape.lastpoint.y - shape.firstpoint.y,
                      shape.lastpoint.x - shape.firstpoint.x)
  radian = radian - (math.pi /2 ) # turn minus 90°
  if (radian > 0):
     degrees = 360 - ( radian  *  360) / ( 2 * math.pi  )
  else:
     degrees = 360 - ((2* math.pi + radian  ) * 360) / ( 2 * math.pi  )
  return degrees
 
# angle =
GetGeographicalDegrees( !SHAPE! )

Now both line features have an angle field and after the spatial join (join one to many – within a distance of XXm) I created a new field that has the difference between angles (AngleDif):

abs(!AngleArcs! - !AngleSegments!)

Next using the python console I used code I found as an answer to a simular problem, which needed to find for every field A with the same value the minimum value of field B:

import arcpy
fc = "C:\\Dropbox\\NewThesis\\databases\\Final.gdb\\ArcosSolSombra_VGA2M_Segments"
data = [row for row in arcpy.da.SearchCursor(fc, ["TARGET_FID", "AngleDif"])]
import collections
ddict = collections.defaultdict(list)
for x,y in data:
   ddict[x].append(y)
minvals = dict((key, min(ddict[key])) for key in ddict.keys())
with arcpy.da.UpdateCursor(fc, ["TARGET_FID", "AngleMin"]) as rows:
    for row in rows:
        row[1] = minvals[row[0]]

The problem I had here is after getting the minvals list the updatecursor is not writing anything in the AngleMin field. 

To finish I'll create and calculate a field to identify when AngleMin is equal to AngleDif.

I know this solution isn't very elegant but as I mentioned before I didn't practice enough after I learnt a bit of Python and now I barely remember anything.

ClintonDow1
Occasional Contributor II

Are you just missing the 'rows.updateRow(row)' in your for loop in the last snippet?

JúlioLuta
New Contributor III

Edit: It worked.

At first I thought it didn't but after I reloaded the file the AngleMin column had the correct values.

Thanks!

0 Kudos