1

I am trying to extract country borders between two points, using Python. The border is a Shapely Polygon.exterior or LinearRing, the start and end points are Shapely Points.

My idea is the following:

  1. Get the closest points or projections for both points along the border

  2. These points split the border into two splines (or LineStrings). Select the shorter one, based on length (ideally reprojected to meters, but it might be overkill).

  3. Decide if the first and/or the last point of the selected spline is needed, to avoid zig-zag like splines at the ends. The given start and end points must be included in the final spline.

  4. Calculated border is: start_point + calculated_spline + end_point


So far I'm a bit stuck with the following concepts:

  • Shapely doesn't support projection on LinearRing, only LineString. How can I calculate the closest points efficiently on a LinearRing? I believe projection isn't needed, as I'll be selecting whether to keep the end points in step 3. Should I just loop over all points of the LinearRing and select the closest one?

  • How can I split a LinearRing into two LineStrings, based on two of their points?

  • Is it true that when comparing relative segment lengths of a LinearRing, I can just measure their length in WGS84?

  • When I have the LineString, should I use .project for deciding weather the endpoints should be included, to avoid zig-zag like splines?

3

I have a solution that works. A LinearRing is also a LineString so use directly a LineString and you can project your points (using Shapely:Coordinate of the closest point on a line)

line = LineString(list(polygon.exterior.coords))
pt1  = line.interpolate(line.project(point1))
pt2  = line.interpolate(line.project(point1))

enter image description here

Now we need to know on what segments of the line are the points and construct the resulting lines from the first point to the second following the original polygon contour (using Determine if shapely point is within a linestring/multilinestring)

def pair(list):
   '''Iterate over pairs in a list -> pair of points '''
   for i in range(1, len(list)):
       yield list[i-1], list[i]

def extract(line, pp1,pp2):
    sol = []
    for seg_start, seg_end in pair(line.coords):
       a = Point(seg_start[0], seg_start[1])
       b = Point(seg_end[0], seg_end[1])
       # after first point found use original segments until second point was reached  
       if len(sol) > 0:  # first point found
            if  LineString([a, b]).distance(pp2) < 1e-8: #second point found
                sol.append(list(LineString([a, pp2]).coords))
                return(MultiLineString(sol))
            else:
               sol.append(list(LineString([a, b]).coords))
       # find first point and construct first line
       if  LineString([a, b]).distance(pp1) < 1e-8:
            sol = [list(LineString([pp1, b]).coords)]  #c.append(list(LineString([pp1, b]).coords))

result = extract(line, pt1,pt2)

enter image description here

enter image description here

But here are some problems:

  • you need to know if the polygon is oriented counter-clockwise or not (in a shapefile, the geometries are counter-clockwise)
  • you need to know the positions of the points in relation to one another (what is the first point ?)
  • When casting from a LinearRing to a LineString it's slightly quicker to do: line = LineString(polygon.exterior). This clones the coordinate sequence using GEOS rather than passing the coordinates via Python. – Snorfalorpagus Jun 27 '15 at 13:44

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.