0

EDIT: Git Repo for sample files https://github.com/tpubben/lineIntersect

I am trying to calculate the line intersection points in x,y coordinates based on a set of intersecting lines crossing one continuous line made up of multiple segments.

The continuous line is represented by a list of tuples as follows where each segment starts with the x/y coordinate of the endpoint of the previous segment:

lineA = [((x1, y1),(x2,y2)), ((x2,y2),(x3,y3))....]

The crossing lines are represented in the same manner, however each is a discrete line (no shared points):

lineB = [((x1, y1),(x2,y2))...]

I am trying to iterate through the continuous line (lineA) and check to see which crossing lines intersect with which segments of lineA.

An example image of what the line intersections would look like is here: Crossing lines

so far I have tried the following:

from __future__ import print_function

def newSurveys(nintyin, injectorin):
    # pull data out of pre-prepared CSV files
    fh = open(nintyin)
    fho = open(injectorin)
    rlines = fho.readlines()
    rlines90 = fh.readlines()

    segA = []
    segB = []
    segA90 = []
    segB90 = []

    for item in rlines:
        if not item.startswith('M'):
            item = item.split(',')
            segA.append((float(item[4]),float(item[5])))#easting northing
            segB.append((float(item[4]),float(item[5])))#easting northing

    segB.pop(0)
    z = len(segA)-1
    segA.pop(z)

    for item in rlines90:
        if not item.startswith('N'):
            item = item.split(',')
            segA90.append((float(item[1]),float(item[0])))#easting northing
            segB90.append((float(item[3]),float(item[2])))#easting northing

    activeWellSegs = []
    injector90Segs = []

    for a, b in zip(segA, segB):
        activeWellSegs.append((a,b))

    for c, d in zip(segA90, segB90):
        injector90Segs.append((c,d))


    if len(activeWellSegs) >= len(injector90Segs):
        lineA = activeWellSegs
        lineB = injector90Segs
    else:
        lineA = injector90Segs
        lineB = activeWellSegs


    for l1 in lineA:        
        for l2 in lineB:
            ##### Use differential equation to calculate line intersections,
            ##### taken from another user's post     
            def line_intersection(line1, line2):
                xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
                ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1])

                def det(a, b):
                    return a[0] * b[1] - a[1] * b[0]

                div = det(xdiff, ydiff)
                if div == 0:
                   raise Exception('lines do not intersect')

                d = (det(*line1), det(*line2))
                x = det(d, xdiff) / div
                y = det(d, ydiff) / div
                return x, y

            print (line_intersection(l1, l2), file=lprint)



    newSurveys('producer90.csv', 'injector.csv') 
haplessgeo
  • 23
  • 4
  • 3
    You have tried that and ... ? What is the actual question? I have little motivation to generate test data, run your code, and then *guess* what your question is. What is the problem with your code? Some actual test data with expected output would be nice. – John Coleman Aug 22 '16 at 19:24
  • Also, do you really want the definition of `line_intersection` to be inside the loop like that? – John Coleman Aug 22 '16 at 19:41
  • I'm pretty new to this...my employer has tossed me into the fire so-to-speak because I have some experience with HTML and CSS (apparently mark-up languages == programming). I don't know if I want that in a nested loop, performance-wise it is probably sub-optimal but raw performance probably doesn't matter too much. How can I attach the test data files for you that is acceptable? – haplessgeo Aug 22 '16 at 21:52
  • Are the number of segments in the `A` list equal to the number of segments in the `B` list, and, if so, does the first segment in `A` intersect the first segment in `B`, etc.? If so, you should be able to replace the nested loop by a single loop. As far as posting data go -- it makes more sense to post examples of the *lists* `lineA` and `lineB`. The file input-output part of your code is a distraction. Separate (into different functions) the logic of extracting the data from the files from the logic of finding the intersection points. – John Coleman Aug 23 '16 at 00:36
  • Unfortunately the data types will never have equal numbers of segments and sometimes there will be more crossing lines than continuous line segments while other times the situation will be reversed. I added a git repo link at the top for anyone who is interested in a sample of the data (and as such the lists). – haplessgeo Aug 23 '16 at 00:56

1 Answers1

-1

Your code looks like it's dealing with a specific set of data (i.e. I have no idea what the "startsWith('N')" and the like are referring to) so I can only answer this question in the abstract.

Try splitting the code into multiple functions that do one specific task, rather than one big function that tries to do everything. You will find it much easier to work with and troubleshoot.

def getScalar(lineSegment):
    return (lineSegment[1][0] - lineSegment[0][0],
            lineSegment[1][1] - lineSegment[0][1])

def doTheyIntersect(lineA, lineB):
    scalarA = getScalar(lineA)
    scalarB = getScalar(lineB)

    s = (-1.0 * scalarA[1] * (lineA[0][0] - lineB[0][0]) + scalarA[0] * (lineA[0][1] - lineB[0][1])) / (-1.0 * scalarB[0] * scalarA[1] + scalarA[0] * scalarB[1])
    t = (scalarB[0] * (lineA[0][1] - lineB[0][1]) - scalarB[1] * (lineA[0][0] - lineB[0][0])) / (-1.0 * scalarB[0] * scalarA[1] + scalarA[0] * scalarB[1])

    if 0.0 <= s <= 1.0 and 0.0 <= t <= 1.0:
        return True
    else:
        return False

lineA = [(x, y), (x1, y1), ...]
lineB = [(x, y), (x1, y1), ...]

for index, segment in enumerate(lineA):
    if index + 1 < len(lineA):
        for index2 in range(0, len(lineB), 2):
             if doTheyIntersect((lineA[index], lineA[index + 1]), (lineB[index2], lineB[index2+1])):
                print("lineB ({0}, {1}) intersects lineA at ({2}, {3})".format(str(lineB[index2]),str(lineB[index2+1]), str(lineA[index]), str(lineA[index + 1]))

This this is the general idea. I got the geometry formulas from:

How do you detect where two line segments intersect?

Community
  • 1
  • 1
TripleD
  • 199
  • 1
  • 15
  • Thanks for getting back to me. I'm working my head through what you've posted and I'm running into a roadblock. After fixing the (remarkably few) syntax errors, I get a traceback that tells me that your getScalar function does not work with tuples. FYI, thanks for reminding me to break it down into functions, much easier to read. – haplessgeo Aug 22 '16 at 22:16
  • For the getScalar function, it need to have a tuple of tuples passed in (e.g. ((1,2), (2,3)) ) with each tuple representing a point of the line segment. There were also a few other mistakes I noticed: forgot to put brackets around the numerator for s and t, the first value in t's formula should be scalarB[0], the if statement is using ">=" instead of "<=", and the print statement at the end is using the wrong formating syntax. I will correct these in the original example. – TripleD Aug 22 '16 at 23:44