Country border between two points

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?

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)) 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, seg_start)
b = Point(seg_end, seg_end)
# 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)  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