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?

  • Welcome to GIS:SE @N'ya! Have you looked at the r.cost tool from GRASS (can also be run via QGIS)?
    – Joseph
    Dec 8, 2015 at 10:51
  • Look at this also gis.stackexchange.com/questions/167287/… Dec 8, 2015 at 13:48
  • @ Joseph: First of all thank you! I´m a step further as I managed to get cost surface rasters already. What I tried to describe is that with common tools (e.g.: saga:leastcostpaths; grass:r.drain) it doesn´t seem to be possible to connect several point with only one processing step.
    – N'ya
    Dec 8, 2015 at 15:27
  • @ Gerado: nice! good link, helped me! what I recognized is that it´s endeed possible to get a cost surface from one point to multiple points (using r.cost). However trying to connect several points within all of them is not possible (as I tried to describe with the hills and valley story)
    – N'ya
    Dec 8, 2015 at 15:30

2 Answers 2


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():
    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
        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')
  • Thank you very much for providing your script and also your previous posts on the topic of LCP! I wasn´t able to check your script yet but it´s next that I´ll do and I´ll let you know if it works for my purpose.
    – N'ya
    Dec 15, 2015 at 10:38
  • @ Anna: just a short but important question about your script. At line 6 "gscript.run_command('r.walk',..." you write "start_coordinates=point_from". Is it right that you just had one starting point and not several starting points? That´s essential for my purpose.
    – N'ya
    Dec 16, 2015 at 14:53
  • I think it's just one point at that particular line. But the whole script takes multiple points. Just run it and see if it's potentially useful for you.
    – Anna
    Dec 17, 2015 at 3:18

Do you have access to ArcGIS?

If you already know in what order you want the LCP to link your points, then I think you could use an iteration with ArcGIS model builder to run your LCP between each points successively:

From A to B
From B to C 

I can provide more information on how to perform it if this solution suits you.

  • 1
    Even though your answer is a possible solution that also came up to my mind already, it´s not exactly what I´m looking for rather what I try to avoid somehow (in order to save resources and time) and find a more elegant way. The most favourable solution would create paths between all point within one step. Furthermore I don´t want to use ArcGIS but what you described would be possible with GRASS too (using r.cost + r.drain iteratively). By the way I´m using PyQGIS (incl. advanced processing toolbox)
    – N'ya
    Dec 8, 2015 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.