grass - Least cost path between various points QGIS

Is there an easy way to find the least cost path (LCP) between various points in an open source GIS?

This question has actually already been asked here: How to find a least cost path through several points? https://gis.stackexchange.com/questions/172894/least-cost-path-through-3-points-start-destination-point-between-them However this post aims to figure out if there´s really no direct way to get favoured results as there were no satisfying answers yet.

There are modules and plug-ins (e.g.: saga:leastcostpaths; grass:r.drain) for least cost path analysis (LCPA) that work well for the LCP between two destination points, but all of them I know and tried so far are not able to compute the LCP between various point.

The reason is that these LCPA-modules compute the path looking for the neighbor pixel with the lowest value of an accumulated cost raster. This works fine as long as there is only one direction (from one point to one other point) to go, however applying these LCPA-modules inputing a cost raster based on various points leads to unsatisfying results as the path got stuck immediatly or after a small amount of pixels.

A cost raster created by various points is like a hilly landscape with ups and down, though LCPA-modules are only able to find the way down they are not able to find the LCP between points, as these point are either located in valleys (so there´s a hills in between) or (at inverted raster) on top of hills (so there are valleys in between where LCP got stuck).

So, I´m looking for a solution how to connect various points by the least cost path based on a cost raster. Does anybody of you has experiences with this topic? Is there an easy way (module/plug-in that already is able to do that job)?

Using the LCPA-modules one cost layer (+LCP) for the distance between each combination of points would be needed. (resources and time consuming)

The created LCP shall build a network suitable for further network analysis.

Edit: In order to specify my question it should be added that LCPA with SAGA GIS (and GRASS GIS probably too) works fine to connect one point with several others but I´m looking for a way to connect several point within themselves. The problem is that GIS tools for LCPA (at least th ones of SAGA and GRASS) are only able to way "downwards" at a cost raster from "peak" where there is the starting point the cost layer was created from, however creating a cost layer with various starting points results in a "hilly" layer where LCPA tools of SAGA and GRASS are not able to find favored paths as they are not able to "climb" but rather got stuck in local minimum pixel values.

Is there any GIS tool able to compute paths between point or polygons based on a cost raster where the points or polygons are located at several local pixel value peaks?

Answer

Here is a script for GRASS GIS 7 I wrote for similar purpose. It takes as input map of vector points and uses r.walk to create cumulative cost surface and creates vector line (least cost path) between each of the points and then computes TSP to get optimal connection between the points using v.net.salesman. You can switch r.walk for r.cost and adjust the parameters based on that but it should work the same.

import grass.script as gscript
import itertools

def trail(elevation, friction, walk_coeff, lambda_, slope_factor, point_from, points_to, vect_paths):
    """Computes least cost path to multiple locations based on walking time"""
    gscript.run_command('r.walk', flags='k', elevation=elevation, friction=friction, lambda_=lambda_, walk_coeff=walk_coeff, slope_factor=slope_factor, start_coordinates=point_from, stop_coordinates=points_to, output='tmp_walk', outdir='tmp_walk_dir')
    for i in range(len(points_to)):
        gscript.run_command('r.drain', flags='d', input='tmp_walk', direction='tmp_walk_dir', output='tmp_drain', drain=vect_paths[i], start_coordinates=points_to[i], overwrite=True)
    gscript.run_command('g.remove', type=['raster', 'vector'], name=['tmp_walk', 'tmp_walk_dir', 'tmp_drain'], flags='f')



def trails_combinations(elevation, friction, walk_coeff, lambda_, slope_factor, points, vector_routes):
    coordinates = gscript.read_command('v.out.ascii', input=points, format='point', separator=',').strip()
    coords_list = []
    for coords in coordinates.split():
        coords_list.append(coords.split(',')[:2])
    combinations = itertools.combinations(coords_list, 2)
    combinations = [list(group) for k, group in itertools.groupby(combinations, key=lambda x: x[0])]
    i = k = 0

    vector_routes_list = []
    for points in combinations:
        i += 1
        point_from = ','.join(points[0][0])
        points_to = [','.join(pair[1]) for pair in points]
        vector_routes_list_drain = []
        for each in points_to:
            vector_routes_list_drain.append('route_path_' + str(k))
            k += 1
        vector_routes_list.extend(vector_routes_list_drain)

        trail(elevation, friction, walk_coeff, lambda_, slope_factor, point_from, points_to, vector_routes_list_drain)
    gscript.run_command('v.patch', input=vector_routes_list, output=vector_routes)


def trails_salesman(trails, points, output):
    gscript.run_command('v.net', input=trails, points=points, output='tmp_net', operation='connect', threshold=10)
    cats = gscript.read_command('v.category', input='tmp_net', layer=2, option='print').strip().split()
    gscript.run_command('v.net.salesman', input='tmp_net', output=output, center_cats=','.join(cats), arc_layer=1, node_layer=2)
    # remove temporary map
    gscript.run_command('g.remove', type='vector', name='tmp_net', flags='f')


if __name__ == '__main__':
    trails_combinations('dem', friction='friction', walk_coeff=[0.72, 6, 0, 6], lambda_=0.5, slope_factor=0, points='markers', vector_routes='route_net')
    trails_salesman(trails='route_net', points='markers', output='route_salesman')