GRASS_movement.qmd

---
title: "Modeling Movement in GRASS"
author: "Michael Barton & Anna Petrasova"
date: 2025-02-10
date-modified: today
lightbox: true
engine: knitter
image: img_movement/thumbnail.webp
format: 
  html:
    embed-resources: true
    toc: true
    code-tools: true
    code-copy: true
    code-fold: false
page-layout: article
categories: [intermediate, advanced, GUI, raster, cost surface, least cost path]
description: Generating a cumulative cost surface and least cost path with *r.walk* and *r.path* to model movement by walking across a landscape.
execute: 
  eval: false
copyright: 
  holder: Michael Barton & Anna Petrasova
  year: 2025
funding: "Creation of this tutorial was supported in part by US National Science Foundation grant FAIN 2303651."

---

GRASS has sophisticated tools to model movement across terrain, including [r.cost](https://grass.osgeo.org/grass-stable/manuals/r.cost.html), [r.walk](https://grass.osgeo.org/grass-stable/manuals/r.walk.html), [r.drain](https://grass.osgeo.org/grass-stable/manuals/r.drain.html), and [r.path](https://grass.osgeo.org/grass-stable/manuals/r.path.html). In this tutorial, we will use ***r.walk*** and ***r.path*** to determine the best walking path to reach a destination, like a hospital.  

::: {.callout-note title="Dataset"}
This tutorial uses one of the standard GRASS sample data sets: flagstaff_arizona_usa. We will refer to place names in that data set, but it can be completed with any of the [standard sample data sets](https://grass.osgeo.org/download/data/) for any region--for example, the [North Carolina data set](https://grass.osgeo.org/sampledata/north_carolina/nc_basic_spm_grass7.zip). We will use the Flagstaff *elevation* DEM, the *hospitals* vector points file (any other vector point file will serve), and the *landuse* raster map.  

This tutorial is designed so that you can complete it using the **GRASS GUI**, GRASS commands from the **console or terminal**, or using GRASS commands in a **Jupyter Notebook** environment.  
:::

::: {.callout-note title="Don't know how to get started?"}
If you are not sure how to get started with GRASS using its graphical user interface or using Python, checkout the tutorials [Get started with GRASS GUI](../get_started/fast_track.qmd) and [Get started with GRASS & Python in Jupyter Notebooks](../get_started/fast_track_grass_and_python.qmd).
:::

## What is a cost surface?

::::: grid
::: g-col-6

-   A **cost surface** is a raster in which each cell represents the “cost” or difficulty to move across landscape.  

-   A **cumulative cost surface** shows the total accumulated cost of moving from a starting point to a location. Cumulative cost surfaces are also used to find a **least cost path** between a location and the starting point.  

-   GRASS [r.walk](https://grass.osgeo.org/grass-stable/manuals/r.walk.html) tool generates a cumulative cost surface using Naismith's rule for walking times where each cell has the value in **seconds** of the time it takes to walk from the start point to that cell.  
:::

::: g-col-6
![cost surface](img_movement/GRASS_movement_0.webp){fig-align="right" width="100%"}  
:::
:::::

# Modeling movement with a cumulative cost surface

## Overview

::::: grid
::: g-col-6

1.  Start with a **DEM of elevation** for topography to determine movement costs.

2.  Create a **friction map** with a value of 0 (or other values for additional movement costs).

3.  Select **start point(s)**.

4.  Create a **cost surface**.

5.  A **least cost** path can then be calculated between any point on the cost surface and the start point.
:::

::: g-col-6
![DEM](img_movement/GRASS_movement_1.webp)
:::
:::::

## Friction map

-   The GRASS [map calculator](https://grass.osgeo.org/grass-stable/manuals/r.mapcalc.html) lets you create a map, alter a map, or combine maps using **map algebra**. The map calculator can be used from the GUI or as a command.  

-   Before you make a friction map, make sure the [**computational region**](https://grass.osgeo.org/grass-stable/manuals/g.region.html) matches the elevation map.  

-   Create the *friction0* map with 0 in every cell.  

:::::: {.panel-tabset group="language"}

#### GUI

::::: grid
::: g-col-6

1.  Add the *elevation* raster to the Layer Manager.  

2.  Right click *elevation* in Layer Manager.  

3.  Select **Set computational region from selected map**.  

4.  Open the **Raster Map Calculator** from the top toolbar ("abacus" icon).  

5.  Fill in the output map name as *friction0* and enter 0 in the expression field.  

6.  Press Run.  
:::

::: g-col-6
![Raster map calculator](img_movement/GRASS_movement_2.webp)
:::
:::::

#### Command line

1.  Set the region to match the elevation map.  

2.  Create a friction map with a value of 0 in all cells that matches the extent and resolution of the elevation map.  

```{bash}
g.region raster=elevation
r.mapcalc "friction0 = 0"
```

#### Python

1.  Set the region to match the elevation map.  

2.  Create a friction map with a value of 0 in all cells that matches the extent and resolution of the elevation map.  

```{python}
gs.run_command("g.region", raster="elevation")
gs.mapcalc("friction0 = 0")
```
::::::

## Choose a start point

Now that we have a DEM map for terrain and a friction map, all we need for a cumulative cost surface is a start point: the point from which to calculate movement costs.  

-   A start point can be a vector point, a raster cell, or even a pair of coordinates.  

-   A cost surface in *r.walk* can have multiple start points.  

-   We’ll make a start point from the Flagstaff Medical Center, found in the *hospitals* vector points file.  

:::::: {.panel-tabset group="language"}

#### GUI

::::: grid
::: g-col-6
Use the **Attribute Table Manager** for *hospitals*:  

1.  Display *hospitals* map by adding it to the Layer Manager from Data Catalog.  

2.  Right click and open the Attribute Table Manager.  

3.  Select "Flagstaff Medical Center" record.  

4.  Right click and select **Extract selected features** from the context menu.  

5.  Name the new map *FMC*.  
:::

::: g-col-6
![data table tool for hospitals vector points map](img_movement/GRASS_movement_3.webp)

![table row context menu](img_movement/GRASS_movement_4.webp){width="60%"}  
:::
:::::

#### Command line

Use the v.extract tool to create a vector point map named *FMC* from the *hospitals* map.  

```{bash}
v.extract input=hospitals type=point where="FACILITY_N = 'FLAGSTAFF MEDICAL CENTER'" output=FMC
```

#### Python

Use the v.extract tool to create a vector point map named *FMC* from the *hospitals* map.  

```{python}
gs.run_command("v.extract", 
                input="hospitals", 
                type="point", 
                where="FACILITY_N = 'FLAGSTAFF MEDICAL CENTER'", 
                output="FMC")
```
::::::

-   Use the *d.vect* tool by double clicking FMC in the layer manager to display the point with a color and size to see it better

![](img_movement/GRASS_movement_5.webp)

## Generate the cumulative cost surface {#generate-cumulative-cost-surface}

::::::::::::: {.panel-tabset group="language"}

#### GUI

Use the *r.walk* tool from the **Raster/Terrain Analysis** menu.  

::::: grid
::: g-col-6

1.  Enter the **elevation map** ( *elevation* ), **friction map** ( *friction0* ), and **name of the cost surface** to create (FMC_cost_seconds).  

![](img_movement/GRASS_movement_6.webp){fig-align="left" width="100%"}  
:::

::: g-col-6
2.  Enter the name of a **directions map** ( *FMC_directions* ) to use for creating a least cost path.  

![](img_movement/GRASS_movement_7.webp){fig-align="left" width="100%"}  
:::
:::::

::::: grid
::: g-col-6
3.  Enter the **start point** ( *FMC* ).  

![](img_movement/GRASS_movement_8.webp){fig-align="left" width="100%"}  
:::

::: g-col-6
4.  Optional: control the **spatial extent** of cost surface.  

![](img_movement/GRASS_movement_9.webp)  
:::
:::::

::::: grid
::: g-col-6
5.  Optional: adjust **Movement parameter** settings.  

![](img_movement/GRASS_movement_10.webp){fig-align="left" width="100%"}  
:::

::: g-col-6
6.  Recommended: select **knight's move** for calculating cost and direction.  

![](img_movement/GRASS_movement_11.webp){fig-align="left" width="100%"}  
:::
:::::

::: {.callout-note title="Tip"}
Click the "copy" button to copy the GRASS command. You can save it in a text file for later reuse or to document your work.  
:::

#### Command line

-   Use the r.walk command to generate the cumulative cost surface.  

```{bash}
r.walk elevation=elevation friction=friction0 output=FMC_cost_seconds outdir=FMC_directions start_points=FMC -k
```

#### Python

-   Use the r.walk command to generate the cumulative cost surface.  

```{python}
gs.run_command("r.walk",
               elevation="elevation",
               friction="friction0",
               output="FMC_cost_seconds",
               outdir="FMC_directions",
               start_points="FMC",
               flags="k")
```
:::::::::::::

## Cumulative cost surface map

Each raster cell value in the cost surface is the time in seconds to walk from FMC to that cell over the terrain of the Flagstaff DEM.  

![](img_movement/GRASS_movement_12.webp)

::: {.callout-note title="Tip"}
Tip: You can display the cumulative cost surface over a shaded relief map in the **layer manager** using the [d.shade](https://grass.osgeo.org/grass-stable/manuals/d.shade.html) tool and a relief map of *elevation* made with [r.relief](https://grass.osgeo.org/grass-stable/manuals/d.shade.html).  
:::

## Movement across a cumulative cost surface {#movement-across-cost-surface}  

You can transform the walking time in seconds to hours by dividing the map by 3600 using the map calculator.  

:::::: {.panel-tabset group="language"}

#### GUI

::::: grid
::: g-col-6

1.  Open the map calculator.  

2.  Enter the name of the new map of walking time in hours: *FMC_cost_hours*.  

3.  Enter FMC_cost_seconds / 3600 into the expressions field.  

4.  Press Run.  
:::

::: g-col-6
![](img_movement/GRASS_movement_13.webp){fig-align="left" width="100%"}  
:::
:::::

#### Command line

```{bash}
r.mapcalc "FMC_cost_hours = FMC_cost_seconds / 3600"
```

#### Python

```{python}
gs.mapcalc("FMC_cost_hours = FMC_cost_seconds / 3600")
```
::::::

------------------------------------------------------------------------

::::: grid
::: g-col-6

-   We can query or filter this hourly cumulative cost surface to areas of equivalent walking time.  

-   For example, in the **layer manager** we can used the *d.rast* tool to show the area within a 2 hour walk of FMC and then set the opacity of the cost surface to 50% to see the underlying terrain.  
:::

::: g-col-6
![](img_movement/GRASS_movement_15.webp)  
:::
:::::

![Shaded area represents terrain reached within a 2 hour walk from FMC](img_movement/GRASS_movement_14.webp)  

# Least cost paths {#least-cost-path}

## Overview

-   We can also plot a **least cost path** (LCP), which is the least costly (shortest time) route between any point on the cumulative cost surface and FMC.  

-   Imagine a stranded hiker NE of Flagstaff who has to walk to FMC.   

-   What path would take the least time?  

![](img_movement/GRASS_movement_16.webp)  

## Generating a least cost path

To create a LCP in GRASS, we will use [r.path](https://grass.osgeo.org/grass-stable/manuals/r.path.html) (also under the Raster/Terrain analysis menu).  

:::::::: {.panel-tabset group="language"}

#### GUI

::::: grid
::: g-col-6

1.  Input the **directions map** (*FMC_directions*) we also made when we created the cumulative cost surface.  

![](img_movement/GRASS_movement_17.webp)  
:::

::: g-col-6
2.  Input the **coordinates of the hiker** as the starting point for the LCP. (We could also use a pre-defined vector point).  

![](img_movement/GRASS_movement_18.webp)  
:::
:::::

:::: grid
::: g-col-6
3.  Specify the **name of the output vector path map** (*LCP_cumulative*).  

![](img_movement/GRASS_movement_19.webp)  
:::
::::

#### Command line

```{bash}
r.path input=FMC_directions format=auto vector_path=LCP_cumulative start_coordinates=477476,13914951
```

#### Python

```{python}
gs.run_command("r.path",
               input="FMC_directions",
               format="auto",
               vector_path="LCP_cumulative",
               start_coordinates=[477476, 13914951])
```
::::::::

## Least cost path generated

Here is the LCP result.  

![](img_movement/GRASS_movement_20.webp)

# Adding a friction map to movement costs

## Overview

-   A **friction map** can be used to incorporate other factors than just terrain into creating a cumulative cost surface and modeling movement.  

-   The value of each cell in a friction map raster is the amount of walking time, in seconds/meter, *in addition to* the time needed because of terrain.  

-   We can **reclassify** the *landuse* map to create a friction map of the amount of extra time it would take to walk through different kinds of land cover.  

![Reclassification of *landuse* to show major land cover types](img_movement/GRASS_movement_21.webp)  

## Reclassifying *landuse* to create a friction map

A standard walking speed across flat terrain is about 5 km/hr = 0.72 sec/m.  

We might then estimate that it would take an additional  

-   3 sec/m to walk through dense pinyon-juniper woodland

-   1 sec/m to cross conifer forest

-   2 sec/m to wander through urban Flagstaff

-   5 sec/m to clamber over lava fields

-   10 sec/m to try and cross water

We can create this friction map by reclassifying the *landuse* map using [*r.reclass*](https://grass.osgeo.org/grass-stable/manuals/r.reclass.html) to assign new friction values to the existing land cover categories.  

:::::: {.panel-tabset group="language"}

## GUI

The [r.reclass](https://grass.osgeo.org/grass-stable/manuals/r.reclass.html) tool is found under the Raster/Change category menu.  

::::: grid
::: g-col-6

1.  Enter *landuse* for the raster to be reclassified.  

2.  Enter *friction_reclassified* for the output reclassified map.  

3.  Enter the reclass rules directly in the text box or from a saved text file. Use and \* symbol to represent everything not covered by the specific reclass rules.  

  > 11 90 95 = 10 water  
  > 21 thru 24 = 2 urban  
  > 31 = 5 lava  
  > 41 thru 43 = 1 conifer forest  
  > 52 = 3 pinyon juniper woodland  
  > \* = 0 no friction  
  
:::

::: g-col-6

![](img_movement/GRASS_movement_22.webp)  
:::
:::::

## Command line

```{bash}
r.reclass input=landuse output=friction_landcover rules=- << EOF
11 90 95 = 10 water
21 thru 24 = 2 urban
31 = 5 lava
41 thru 43 = 1 conifer forest
52 = 3 pinyon juniper woodland
* = 0 no friction
EOF
```

## Python

```{python}
rules = """\
11 90 95 = 10 water
21 thru 24 = 2 urban
31 = 5 lava
41 thru 43 = 1 conifer forest
52 = 3 pinyon juniper woodland
* = 0 no friction
"""
gs.write_command("r.reclass",
                 input="landuse",
                 output="friction_landcover",
                 rules="-",
                 stdin=rules)
```
::::::

This creates a new friction map named *friction_landcover*.  

![The reclassified friction map (*friction_landcover*)](img_movement/GRASS_movement_23.webp)  

## Modifying a cumulative cost surface with a friction map

We can now create a new cumulative cost surface using this new friction map and convert it from seconds to hours, as we did before.  

::: {.panel-tabset group="language"}

#### GUI

1. Follow the procedures in the [Generate the cumulative cost surface
section](#generate-cumulative-cost-surface) and substitute the new *friction_landcover* map for the *friction0* map used previously.  

2. Convert the cumulative cost surface to hours instead of seconds following the procedures in the [Movement across a cumulative cost surface section]({#movement-across-cost-surface}).  

#### Command line

```{bash}
r.walk elevation=elevation friction=friction_landcover output=FMC_vegcost_seconds outdir=FMC_vegcost_directions start_points=FMC -k
r.mapcalc "FMC_vegcost_hours = FMC_vegcost_seconds / 3600"
```

#### Python

```{python}
gs.run_command("r.walk",
               elevation="elevation",
               friction="friction_landcover",
               output="FMC_vegcost_seconds",
               outdir="FMC_vegcost_directions",
               start_points="FMC",
               flags="k")
gs.mapcalc("FMC_vegcost_hours = FMC_vegcost_seconds / 3600")
```
:::

![cumulative cost surface with additional friction from land cover](img_movement/GRASS_movement_24.webp)  

## Least cost paths and friction

We can create a new LCP from the stranded hiker to FMC over terrain where land cover also affects the cost of movement.  

::: {.panel-tabset group="language"}

#### GUI

1. Follow the procedures in the [Least cost path section](#least-cost-path), substituting the new direction map made along with the cumulative cost surface with land cover friction map.  

2. Give this new LCP the name *LCP_vegcost* to distinguish from the previous LCP.  

#### Command line

```{bash}
r.path input=FMC_vegcost_directions format=auto vector_path=LCP_vegcost start_coordinates=477476,13914951
```

#### Python

```{python}
gs.run_command("r.path",
               input="FMC_vegcost_directions",
               format="auto",
               vector_path="LCP_vegcost",
               start_coordinates=[477476, 13914951])
```
:::

-   It takes more time to reach FMC if the hiker has to navigate dense vegetation as well as terrain.  

-   In the map below, terrain colored by land cover, the original LCP with terrain only is shown by the blue line. The LCP with a land cover friction map is shown by the heavier yellow line.  

![comparing LCP with terrain only and with land cover friction](img_movement/GRASS_movement_25.webp)