Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Determine the "left" and "right" side of a split shapely geometry

Tags:

python

shapely

My question is: how can I determine which of the Aside and Bside sides of an already split rotated rectangular geometry are the "left" and "right" side of an arbitrary LineString that splits that geometry?

For the purposes of this problem, "left" and "right" are defined as the left-hand and right-hand side of a LineString splitter when "walking" from node to node, in order.

I've created this function for splitting any arbitrary shapely geometry (non collection) into two sides - the "left" and the "right":

import shapely.geometry as geo
import shapely.ops as ops

def splitLR(geom, splitter):
    """Split a geometry into a 'left' and 'right' side using the shapely API"""
    if not isinstance(splitter, geo.LineString):
        raise TypeError("The splitter must be a LineString")
    if not splitter.is_simple:
        raise ValueError("Only simple splitter objects allowed")
    if hasattr(geom, "__iter__"):
        raise ValueError("Geometry collections not allowed")
    geom_extents = geo.GeometryCollection([geom, splitter]).minimum_rotated_rectangle
    sides = ops.split(geom_extents, splitter)
    try:
        Aside, Bside = sides
    except TypeError:
        # only 1 result - rotated rectangle wasn't split
        if len(ops.split(geom,splitter)) == 1:
            # geom isn't split by splitter
            raise ValueError("the splitter does not appear to split the geometry")
        else:
            # splitter too small for algorithm
            raise ValueError("the splitter must extend beyond minimum_rotated_rectangle "
                             "of the combined geometry")
    # determine which is Lside and Rside here
    Lside,Rside = get_LRsides(Aside, Bside, splitter)
    return tuple(side.intersection(geom) for side in (Lside, Rside))

The idea of the above is illustrated in the notebook linked here (same link as above):

http://nbviewer.jupyter.org/urls/dl.dropbox.com/s/ll3mchnx0jwzjnf/determine%20left-right%20split.ipynb

To summarize: sides A and B are the two sides of the minimum_rotated_rectangle surrounding the geom geometry and the splitter line string together. When side.intersection(geom) is executed, the result is the portion of the originally given geometry geom contained in that side.

Notes:

  • it is possible- for oddly shaped "potato" type objects- for this intersection to result in multiple objects on one or both sides (see the nbviewer example)
  • I have created my own function here (rather than using ops.split) because the ops.split function just returns a "bag" of split objects, and there isn't a way to determine which side they are on (that I'm aware of)

Currently my call to get_LRsides just executes this function, which is obviously worthless:

def get_LRsides(Aside, Bside, splitter):
    """Determine the 'left' and 'right' sides of an already split geometry"""
    return Aside,Bside

How can I successfully label A and B as "left" and "right"?

like image 255
Rick supports Monica Avatar asked Sep 05 '25 02:09

Rick supports Monica


1 Answers

This could work:

  1. Form a LinearRing with the end-points of the splitter and a point in Aside
  2. Apply object.is_ccw
  3. If it returns True, Aside is to the left of the splitter.

https://shapely.readthedocs.io/en/stable/manual.html#object.is_ccw

like image 171
yang5 Avatar answered Sep 07 '25 21:09

yang5