Why LineString can't be splitted to 2 parts？

[forestbat]

I wrote this code：

nearest_p = nearest_points(source_coord, origin_line)[1]
# source_coord: POINT (132.775278 45.298611)
# origin_line: LINESTRING (132.785417 45.310417, 132.772917 45.297917, 132.764583 45.297917, 132.760417 45.29375, 132.66875 45.29375)

And nearest_p is POINT (132.7744445 45.29944449999999), but when I run this code:
min_lines = split(origin_line, nearest_p)
length of min_lines is 1 and in fact, nearest_p is not on origin_line.
So what happened?

[brendan-ward]

I'm pretty sure this is due to precision / robustness issues; it is difficult to produce a point that is exactly on the line, so you need to use a different approach to split this. As noted in the docstring for split you may need to use lower precision to help make sure the

Note that you can also interpolate the point into the line to calculate its approximate coordinate:

>>> import shapely
>>> shapely.line_interpolate_point(origin_line, shapely.line_locate_point(origin_line, source_coord))
<POINT (132.774 45.299)>

From there, you can use an algorithm that is similar to that used by split but bypass the first check that the point intersects the line (which is returning False in this case).

[forestbat]

@brendan-ward Thanks for very much，but the solution can't solve my problem completely.
I change my code like this：

def sure_nearest_point(source_coord: Point, origin_line: LineString):
    nearest_p = nearest_points(source_coord, origin_line)[1]
    distance_on_line = origin_line.project(source_coord)
    tolerance = 1e-5 
    point_src = origin_line.coords[0]
    point_dst = origin_line.coords[-1]
    if (((nearest_p.x == point_src[0]) & (nearest_p.y == point_src[1])) |
       ((nearest_p.x == point_dst[0]) & (nearest_p.y == point_dst[1]))) | (distance_on_line < tolerance)\
            | (distance_on_line > shapely.length(origin_line) - tolerance):
        if len(origin_line.coords) > 2:
            rand_coord = random.randint(0, len(origin_line.coords) - 1)
            nearest_point = Point(origin_line.coords[rand_coord])
        else:
            nearest_point = Point((point_src[0] + point_dst[0]) / 2, (point_src[1] + point_dst[1]) / 2)
    else:
        nearest_point = shapely.line_interpolate_point(origin_line, shapely.line_locate_point(origin_line, source_coord))
    return nearest_point

And there is no other change in my code，but when I use other points' shapefile which is different from points.zip in this discussion，the program still crashed and adjusting tolerance has no effect.

[brendan-ward]

The tolerance parameter I suggested is for determining if you should split the line after finding the nearest point, not for determining the nearest point.

[forestbat]

I wrote a test for ensuring tolerance：

def test_tolerance():
    gpd_min_line_df = gpd.sjoin_nearest(gpd_df_node, gpd_df_network, 'left')
    gpd_min_line_df = gpd_min_line_df.loc[~gpd_min_line_df.index.duplicated()]
    index_right = gpd_min_line_df['index_right'].to_numpy()
    origin_geom_array = gpd_df_network.geometry.to_numpy()
    project_list = []
    for index in range(0, len(gpd_df_node)):
        point = gpd_df_node.geometry[index]
        right = index_right[index]
        near_point = sure_nearest_point_old(point, origin_geom_array[right])
        project_list.append(origin_geom_array[right].project(near_point))
    print(project_list)

def sure_nearest_point_old(source_coord: Point, origin_line: LineString):
    nearest_p = nearest_points(source_coord, origin_line)[1]
    point_src = origin_line.coords[0]
    point_dst = origin_line.coords[-1]
    if ((nearest_p.x == point_src[0]) & (nearest_p.y == point_src[1])) | \
            (nearest_p.x == point_dst[0]) & (nearest_p.y == point_dst[1]):
        if len(origin_line.coords) > 2:
            rand_coord = random.randint(0, len(origin_line.coords) - 1)
            nearest_point = Point(origin_line.coords[rand_coord])
        else:
            nearest_point = Point((point_src[0] + point_dst[0]) / 2, (point_src[1] + point_dst[1]) / 2)
    else:
        nearest_point = shapely.line_interpolate_point(origin_line, shapely.line_locate_point(origin_line, source_coord))
    return nearest_point

And project_list looks like this：

[0.04432260408057233, 0.04218472542212115, 0.0637690182941307, 0.004222999999996092, 0.015517458313150477, 0.09440157625101754, 0.0003060000000019158, 0.0772293487644716, 0.05003792933941874, 0.025677669529667435, 0.032753191317823255, 0.08558036702676952, 0.0105829999999969, 0.10315567835199253, 0.06703041916308822, 0.019416231165248918, 0.006528000000002976, 0.06037599046457332, 0.018421027914405343, 0.010027000000000898, 0.02226096246530259, 0.05426673554460159, 0.047198147565396996, 0.03337196246530016, 0.06307111799533686, ……]

You can realize 0.04 degree is too big to "ignore" them in a tolerance, so what's the problem and how to solve?

[brendan-ward]

I don't understand why you are doing this block:

    if ((nearest_p.x == point_src[0]) & (nearest_p.y == point_src[1])) | \
            (nearest_p.x == point_dst[0]) & (nearest_p.y == point_dst[1]):
        if len(origin_line.coords) > 2:
            rand_coord = random.randint(0, len(origin_line.coords) - 1)
            nearest_point = Point(origin_line.coords[rand_coord])
        else:
            nearest_point = Point((point_src[0] + point_dst[0]) / 2, (point_src[1] + point_dst[1]) / 2)

which, if I'm interpreting it correctly, is basically: if the nearest point is the first or last point of the line and the line has more than 2 points, pick a random vertex on that line, and if the line only has 2 points, use the midpoint on the line

Both of which offset your nearest point from the actual nearest point. If the nearest point on the line to your original point is an endpoint (first or last point) on the target line, there is no reason to split the line: the line is already split at that point!

Instead, I think you could replace the call to sure_nearest_point with

near_point = shapely.line_interpolate_point(
    origin_geom_array[right], shapely.line_locate_point(origin_geom_array[right], point)
)

Because that also finds the nearest coordinate.

This might make it a little more clear:

point = Point(0,0)
line = LineString([(1,1), (2,2)])

shapely.line_interpolate_point(line, shapely.line_locate_point(line, point))
# <POINT (1 1)>

point = Point(1.5,1.25)
shapely.line_interpolate_point(line, shapely.line_locate_point(line, point))
# <POINT (1.375 1.375)>

Using this nearest point, you then need to decide if you are going to split the line or if it is already effectively split because the nearest point to your original point is the endpoint of the line. Since testing equality of floating point numbers is problematic, one way to approach this is to calculate the position on the line and use a tolerance to make the decision:

tolerance = 1e-10

nearest_point = shapely.line_interpolate_point(line, shapely.line_locate_point(line, point))

pos = shapely.line_locate_point(line, nearest_point)
if pos > tolerance and pos < shapely.length(line) - tolerance:
    # split the line...
else:
    # use the original line

This ensures that you only split the line, and add a new coordinate, if that new coordinate is between the endpoints of the line.

If you wanted, you could make this more complex by checking to see if the nearest point is within a tolerance of one of the existing vertices on the line other than the endpoint. E.g., if your line has 3 vertices and the point is very very close to the middle vertex, then reuse that existing vertex as the split point instead.

If you are doing this for many lines and points, it would benefit from taking a vectorized approach - but I'd recommend getting a solid handle on the pure python approach first. Some ideas at: pygeos/pygeos#352

Some functionality in this space has been added to GEOS recently but I don't think that has been implemented in Shapely yet.

[forestbat]

Yes now I replaced my code with your method：

def test_tolerance():
    gpd_min_line_df = gpd.sjoin_nearest(gpd_df_node, gpd_df_network, 'left')
    gpd_min_line_df = gpd_min_line_df.loc[~gpd_min_line_df.index.duplicated()]
    index_right = gpd_min_line_df['index_right'].to_numpy()
    origin_geom_array = gpd_df_network.geometry.to_numpy()
    project_list = []
    # Time complexity is O(m), m is amount of points
    for index in range(0, len(gpd_df_node)):
        point = gpd_df_node.geometry[index]
        right = index_right[index]
        near_point = shapely.line_interpolate_point(origin_geom_array[right],
                                                    shapely.line_locate_point(origin_geom_array[right], point))
        project_list.append(origin_geom_array[right].project(near_point))
    print(project_list)

But project list is still like this, too big:

0.04432260408057233, 0.04218472542212115, 0.0637690182941307, 0.004222999999996092, 0.015517458313150477, 0.09440157625101754, 0.0003060000000019158, 0.0772293487644716, 0.05003792933941874, 0.025677669529667435, 0.032753191317823255, 0.08558036702676952, 0.0105829999999969, 0.10315567835199253, 0.06703041916308822, 0.019416231165248918, 0.006528000000002976……

My algorithm to find nearest lines comes from here.

[forestbat]

Now I use a new method to solve this problem：

def line_min_dist(gpd_node_df: GeoDataFrame, gpd_line_df: GeoDataFrame):
    gpd_min_line_df = gpd.sjoin_nearest(gpd_node_df, gpd_line_df, 'left')
    gpd_min_line_df = gpd_min_line_df.loc[~gpd_min_line_df.index.duplicated()]
    index_right = gpd_min_line_df['index_right'].to_numpy()
    origin_geom_array = gpd_line_df.geometry.to_numpy()
    geom_array = np.delete(origin_geom_array, index_right)
    for index in range(0, len(gpd_node_df)):
        point = gpd_node_df.geometry[index]
        right = index_right[index]
        pre_near_point = sure_nearest_point(point, origin_geom_array[right])
        line = LineString([point, pre_near_point])
        if line.length != 0:
            ex_line = get_extrapolated_line(point, pre_near_point, 1 / line.length)
            min_lines = split(origin_geom_array[right], ex_line)
        else:
            min_lines = split(origin_geom_array[right], pre_near_point)
        geom_array = np.append(geom_array, (list(min_lines.geoms))[1])
    return geom_array


def sure_nearest_point(source_coord: Point, origin_line: LineString):
    nearest_p = nearest_points(source_coord, origin_line)[1]
    point_src = origin_line.coords[0]
    point_dst = origin_line.coords[-1]
    point_neigh = origin_line.coords[1]
    if ((nearest_p.x == point_src[0]) & (nearest_p.y == point_src[1])) | ((nearest_p.x == point_dst[0])
                                                                          & (nearest_p.y == point_dst[1])) | \
            ((source_coord.x == nearest_p.x) & (source_coord.y == nearest_p.y)):
        nearest_point = Point((point_src[0] + point_neigh[0]) / 2, (point_src[1] + point_neigh[1]) / 2)
    else:
        nearest_point = shapely.line_interpolate_point(origin_line,
                                                       shapely.line_locate_point(origin_line, source_coord))
    return nearest_point


def get_extrapolated_line(source_coord, coord, extrapolate_ratio):
    source = Point(source_coord)
    p2 = Point(coord)
    c = (source.x + extrapolate_ratio * (p2.x - source.x), source.y + extrapolate_ratio * (p2.y - source.y))
    return LineString([source, c])