1

Is there an easy solution for recording events for a line crossing borders of zones in a polygon layer?

I need to implement this in python using QGIS, thus pyqgis. The image attached illustrates what I am looking for. The red line is a linestring geometry with start and end point. Whenever the line enters a zone (i.e. crosses the polygon border) an event is recorded being the reference of the zone (A to D) and the distance from the start of the line, measured along the linestring's geometry (distances in the image are estimates only). I can get by the first "event", where the line already starts inside zone A. But is there an easy way to get the other events? I could cut up the line in segments using the zones, but my goal is to record events to prepare for linear referencing, so I try to avoid cutting up the line.

enter image description here

1

I was able to figure this out myself using the Shapely library, like so:

from shapely.geometry import LineString, Point, Polygon

epsilon = 0.001
def is_between(a, b, c):
    crossproduct = (c.y - a.y) * (b.x - a.x) - (c.x - a.x) * (b.y - a.y)
    if abs(crossproduct) > epsilon : return False # (or != 0 if using integers)
    dotproduct = (c.x - a.x) * (b.x - a.x) + (c.y - a.y) * (b.y - a.y)
    if dotproduct < 0 : return False
    squaredlengthba = (b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y)
    if dotproduct > squaredlengthba: return False
    return True

def get_lin_ref(linestring, point):
    lin_ref = 0
    for i in range(len(linestring.coords) - 1):
        point0 = Point(linestring.coords[i])
        point1 = Point(linestring.coords[i + 1])
        if is_between(point0, point1, point):
            lin_ref += point0.distance(point)
            return lin_ref
        else:
            lin_ref += point0.distance(point1)

def get_events(linestring, polygon):
    intersect = linestring.intersection(polygon)
    events = []
    if intersect.geom_type == "GeometryCollection":
        for geometry in intersect:
            if geometry.geom_type == "MultiLineString":
                for intersect_part in geometry:
                    events.append((
                            get_lin_ref(linestring, Point(intersect_part.coords[0])),
                            get_lin_ref(linestring, Point(intersect_part.coords[-1]))
                        ))
            elif geometry.geom_type == "LineString":
                events.append((
                        get_lin_ref(linestring, Point(geometry.coords[0])),
                        get_lin_ref(linestring, Point(geometry.coords[-1]))
                    ))
    elif intersect.geom_type == "MultiLineString":
        for intersect_part in intersect:
            events.append((
                    get_lin_ref(linestring, Point(intersect_part.coords[0])),
                    get_lin_ref(linestring, Point(intersect_part.coords[-1]))
                ))
    elif intersect.geom_type == "LineString":
        events.append((
                get_lin_ref(linestring, Point(intersect.coords[0])),
                get_lin_ref(linestring, Point(intersect.coords[-1]))
            ))
    return events

# Example:
line = LineString([(1, 1), (1, 3), (3, 3), (5, 3), (7, 3), (9, 3), (9, 5)])
ext = [(2, 1), (8, 1), (8, 5), (2, 5), (2, 1)]
int = [(4, 2), (6, 2), (6, 4), (4, 4), (4, 2)]
polyA = Polygon(ext, [int])
polyB = Polygon(int)

polyA_events = get_events(line, polyA)
polyA_events = [["polyA"] + list(tup) for tup in polyA_events]

polyB_events = get_events(line, polyB)
polyB_events = [["polyB"] + list(tup) for tup in polyB_events]

all_events = polyA_events + polyB_events
all_events.sort(key=lambda event: event[1])

print all_events
# prints: [['polyA', 3.0, 5.0], ['polyB', 5.0, 7.0], ['polyA', 7.0, 9.0]]

Any improvement would be very much appreciated. Thanks.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.