movecost (R package) vignette

movecost R package Dr Gianmarco Alberti - University of Malta 23/08/2022 Introduction to the package’s vignette This vignette aims at showing the use of the current version of the movecost package and of its functions. To hopefully enhance clarity, it is organised as a sequence of tasks. In-built datasets are used throughout this document so that the users may easily familiarise themselves with the functions and their parameters. For more details about each function, its parameters, for the values returned, for the export options, and for relevant literature, users are referred to the package’s help documentation. Task 1. Given an origin location, how can I calculate a cost surface accumulated around that point? This can be easily achieved with the movecost() function. We need a Digital Terrain Model and the origin location. In this example, we can use some in-build datasets. Assuming that we want to calculate the accumulated cost of moving outward from the origin, framing the cost in terms of walking time, we can run what follows: result <- movecost(dtm=volc, origin=volc.loc, move=16, funct="t", time="m", breaks=2) In the above code, we have defined the cost function to be used (the Tobler’s hiking function in this case) via the funct parameter, the time unit to use (minutes) via the time parameter, and the isolines/isochrones interval (2 minutes) via the breaks parameter. Task 2. How can I calculate a cost surface accumulated around an origin, using a terrain factor to account for a given terrain type? This task can be accomplished using virtually the same command as above, but using the N parameter to express the terrain factor. result <- movecost(dtm=volc, origin=volc.loc, move=16, funct="t", time="m", breaks=2, N=1.19) In the above example, a terrain factor N of 1.19 has been used, corresponding to loose beach sand. It is apparent from the above image that the walking time expands because moving on sand reduces the walking speed (i.e., increase the pace). Task 3. Can cost be expressed in terms of metabolic energy expenditure rather than walking time? Yes. Many energetic cost functions are implemented in the movecost package. The following example is the same as above, but uses a different cost function. result <- movecost(dtm=volc, origin=volc.loc, move=16, funct="pcf", W=60, L=5) Here, the cost is expressed in terms of energy spent while walking according to the Pandolf et al.’s cost function incorporating a correction factor. The parameters W and L define the body weight anf the load carried (both in Kg). Note that the V parameter (walking speed, in m/s) has not been entered and the deafult value 1.2 is internally used. In this case, the walking speed is considered uniform across the entire study area (more on this below). Task 4. When using some energetic cost functions, can the walking speed be internally worked out instead of using a single value for the entire study area? Yes. The following example is similar to the previous one, with one exception. result <- movecost(dtm=volc, origin=volc.loc, move=16, funct="pcf", W=60, L=5, V=0) Here, the parameter V has been set to 0 . In this case V is internally derived using the Tobler’s hiking function, so the walking speed will vary as a function of the terrain slope. Task 5. Can the so-called cognitive slope be used when calculating cost surfaces? Yes. The following example uses the same settings of Task 1 above, but uses the cognitive slope instead. result <- movecost(dtm=volc, origin=volc.loc, move=16, funct="t", time="m", breaks=2, cogn.slp=TRUE) In this case, the parameter cogn.sl has been set to TRUE in order to use the cognitive slope. As a result, following Pingel TJ (2013), Modeling Slope as a Contributor to Route Selection in Mountainous Areas, in Cartography and Geographic Information Science, 37(2), 137-148, before being used for the cost calculation positive slope values are preliminarily multiplied by 1.99, while negative slope values are multiplied by 2.31. Task 6. Can least-cost path(s) between an origin and one or multiple destinations be calculated? Yes. The following example uses the same dataset used so far, plus another in-built dataset ( destin.loc ) that can be readily used. result <- movecost(dtm=volc, origin=volc.loc, destin=destin.loc, move=16, funct="t", time="m", breaks=2, oneplot= FALSE) In this case, we have fed the destination locations into the function via the destin parameter. Two plots are rendered. One (not shown here) representing the accumulated cost surface around the origin; another (the above) representing the least-cost paths (hereafter LCPs). The destinations feature labels indicating the cost at those locations. The cost, in this example, is expressed in terms of walking time on the basis of the Tobler’s function for on-path hiking. The walking-time distance at the destination locations is expressed in sexagesimal format (hh:mm:ss). The two mentioned plots can be arranged in one single visualisation (i.e., side-by-side) or separately by using the parameter oneplot (which in this case has been set to FALSE ). Task 7. Can LCPs back to the origin be calculated? Yes. The following code is the same as for Task 6, but uses the parameter return.base . result <- movecost(dtm=volc, origin=volc.loc, destin=destin.loc, move=16, funct="t", time="m", breaks=2, oneplot= FALSE, return.base=TRUE) As shown in the above figure, dashed lines represent the LCPs from the destinations back to the origin. This is easily achieved setting the parameter return.base to TRUE . Task 8. Can a least-cost corridor be calculated? Yes. This can be easily achieved with the movecorr() function. The following code calculates the least-cost corridor between the origin location used in the preceding examples and one of the destination locations used in task 6 and 7. result <- movecorr(dtm=volc, a=volc.loc, b=destin.loc[2,], funct="wcs") The above image shows the least-cost corridor raster between the two locations; the cost of a cell is the total cost to reach it from the two locations. In the rendered plot, least-cost paths from A to B and viceversa are also shown. In this example, the cost function used (fed via the funct parameter) is the cost function for wheeled vehicle. The critical value of the slope for this function is set by default to 10 (percent); it can be modified using the sl.crit parameter. Note that energetic cost functions, terrain factor , cognitive slope, and barriers can be used in movecorr() as well. Task 9. Can I compare LCPs based on different cost functions? Yes. This can be easily achieved with the movecomp() function. The following code will calculate and plot LCPs based on three different cost functions expressing cost in terms of walking time. To enhance readability, only three destination locations are used. Notice the add.chart parameter, set to TRUE , which will be commented on later on. result <- movecomp(dtm=volc, origin=volc.loc, destin=destin.loc[c(2,4,6),], choice=c("t", "ma", "ug", "gkrs"), mo ve=16, add.chart=TRUE) The cost functions to compare have been fed into the choice parameter using the character vector c("t", "ma", "ug", "gkrs") . In the above image, LCPs are rendered on the DTM. Different line types indicate the cost-function used for each path. The LCPs back to the origin can be rendered setting the parameter return.base to TRUE . Note that energetic cost functions, terrain factor, cognitive slope, and barrier can be used in movecomp() as well. The above-mentioned add.chart parameter will produce two additional plots, one of which is shown below: The above image shows the boxplots portraying the distributon of the length of the LCPs by cost function; it is apparent that one of the used cost functions (namely, gkrs) produces longer LCPs compared to the other three cost functions. Another boxplot (not shown here) rendered by the function shows the distribution of the cost by cost function. Task 10. I would like to use either of the functions described so far, but I do not have an elevation raster; what can I do? All the functions so far described can download elevation data from online sources using (under the hood) functions out of the elevatr package by Jeffrey Hollister and colleagues. If you have a polygon ( SpatialPolygonDataFrame ) representing the study area, movecost can download elevation data for the area enclosed by that polygon. The zoom level of the downloaded elevation data (i.e., its resolution) is controlled by the parameter z , which is set to 9 by default (a trade off between resolution and download time). The following example is similar to the preceding Task 9 in that it compares LCPs based on different cost functions, but aquires elevation data for an area around Mt Etna (Sicily, Italy). result <- movecomp(studyplot=Etna_boundary, origin=Etna_start_location, destin=Etna_end_location, choice=c("t", "wcs"), move=16) The description of the above image is along the same lines of the description of the plot rendered for Task 9. Elevation data have been aquired for the area enclosed by the polygon fed via the studyplot parameter. The datasets used as origin and destination locations are in-built in the movecost package. For more information about the elevation data resolution per zoom level, and to know what is sourced at what zoom level, see the links provided in the movecost ’s help documentation. Task 11. Can I calculate and plot a specific walking-cost boundary around one or multiple locations? For instance, can I calculate a 45-minute walking boundary around some locations? Yes. This can be easily achieved with the movebound() function. Consider the following code: result <- movebound(studyplot = Etna_boundary, origin=Etna_end_location, cont.value=45, time="m", funct="tofp", a dd.geom=TRUE) In the above image, the 45-minute walking boundaries (based on the Tobler’s off-path hiking function) are plotted on the DTM around the three locations. Since the studyplot parameter is used, the function downloads elevation data from online sources. By setting the add.geom parameter to TRUE , the perimeter and area of the calculated boundaries are stored as two new variables appended to a copy of the input origin locations dataset, which is returned by the function. Needless to say, energetic cost functions, terrain factor and cognitive slope can be used in movebound() as well. Task 12. Can I carry out a cost allocation analysis? Yes. This can be easily achieved with the movealloc() function, which provides the facility to carry out cost allocation analysis. Given a number of origin locations, a cost allocation raster is produced; each cell of the cost allocation raster is given an integer indicating to which origin a cell is closer in terms of cost. Needless to say, the cost can be conceptualized in terms of either walking time or energy expenditure, and is function of the terrain slope. Consider the following code: result <- movealloc(dtm=volc, origin=destin.loc[c(3,7,9),], funct="a") In the above image, the cost allocation raster is rendered on the input DTM ( dtm=volc ); only three locations are used for the analysis ( origin=destin.loc[c(3,7,9),] ). The cost function used (fed via the funct parameter) is the Ardigo et al.’s for metabolic cost; the parameters of this cost function are left as they are by default. The study area is divided into three zones; each zone is given a colour to indicate to which origin each cell is closer in terms of walking metabolic cost. For instance, the cell that are given yellow are closer to the origin location of the yellow area than to any other origin location. Note that energetic cost functions, terrain factor and cognitive slope can be used in movealloc() as well. Task 13. In the cost allocation plot, can I also show isolines inside each cost allocation area? Yes. The following code is similar to the one used in the previous task, with one exception: result <- movealloc(dtm=volc, origin=destin.loc[c(3,7,9),], funct="a", isolines=TRUE) In the above image, isolines are plot within each allocation zone. This is easily achieved by setting the parameter isolines to TRUE . The isolines’ (i.e., contours’) labels can be deactivated by setting the cont.lab parameter to FALSE ( TRUE is the default value). Task 14. Can I calculate a network of LCPs (in other words, LCPs between multiple origins)? Yes. It is possible to calculate two different types of networks: one where each origin location is connected to all the others locations; one where only pairs of neighboring locations are connected. In other words, in the latter case, each location is connected to the location that is the nearest in terms of walking cost, either in terms of time or energy (or abstract cost), according to the selected cost function. Either network type is produced by the movenetw() function: result1 <- movenetw(dtm=volc, origin=destin.loc, funct="t", netw.type="allpairs") The above image shows the first type of network, which is produced setting the netw.type parameter to allpairs (which is the default value). LCPs connect every point location (represented by the red dots) to each other location (the on-path Tobler’s hiking function is used). The second type of network is produced with the code below: result2 <- movenetw(dtm=volc, origin=destin.loc, funct="t", netw.type="neigh") text(destin.loc, labels=c(1:length(destin.loc)), cex=0.6, adj=c(-1,-1)) The second line of code adds labels indicating the locations number, and is instrumental to the comments provided later on when reviewing the cost matrix. In this second type of network (note that netw.type parameter has been set to neigh ), every location is connected by a LCP to the location that is the closest in terms of time spent (or energy consumed, or abstract cost incurred) while walking. In this example, since we are using the on-path Tobler’s function, the connected locations are closest in terms of walking-time distance. We can inspect the matrix storing all the pairwise walking-time distances by using what follows: result2$cost.matrix.min From the matrix, first and foremost we can see that the values are not symmetric about the diagonal for the very fact that we are dealing here with anisotropic costs (i.e., the cost is direction-dependent: AtoB is not necessarily equal to BtoA). Secondly, reading the matrix off row-wisely, we can see for example that location 9 is the closest to location 1, and they are therefore connected by a LCP. However, the location that is closest to location 9 is not location 1 but location 7, and they are therefore linked. It is worth noting that location 5 is connected to location 6 and viceversa because they are the closest to one another; this is the reason why between the two locations we can see two partially overlapping LCPs (one from 5 to 6, one from 6 to 5). Finally, if we consider location 3, it is worth noting that location 4 is the closest in terms of planar distance, but location 2 is indeed the closest to location 3 in terms of walking-time distance (5.64 vs. 6.14 minutes). They are therefore linked by a LCP. Task 15. Can I calculate the density of a LCPs network? Yes. The density of the LCPs network can be calculated for the first type of network described earlier on, namely the one connecting all the locations. The following code is similar to the one used for the preceding task, but employs the lcp.dens parameter, which is set to TRUE result <- movenetw(dtm=volc, origin=destin.loc, funct="t", netw.type="allpairs", lcp.dens=TRUE) The above command will produce two plots: one is the LCPs network plot (see Task 14 above); the second is the one shown here. It features an hillshade visualization, on top of which the LCPs density raster is rendered. The density of overlapping LCPs is expressed as percentage. The percentages are calculated in relation to the maximum number of LCPs passing through the same cell. A raster expressing density in terms of counts is not plotted but returned by the function. Task 16. My DTM has irregular margins, i.e. the study area features gulfs and inlets, and the layer contains NoData values corresponding to the sea. The calculated LCP crosses the sea: how can I prevent that? As of version 1.8, the user can use the irregular.dtm parameter to prevent the LCP(s) the cross the sea. Compare the two following examples: resultA <- movecost(malta_dtm_40, springs[5,], springs[15,]) In the above example, the DTM represents an island (namely, Malta), and it features irregular margins with gulfs and inlets. The area corresponding to the sea is given NoData. As apparent from the figure, the LCP does cross the sea. Instead, by using the following code resultB <- movecost(malta_dtm_40, springs[5,], springs[15,], irregular.dtm=TRUE) the returned plot shows a LCP that, at some point, follows the coastline. Internally, what movecost() does is to generate a polygon vector layer from the DTM and to use the polygon as a mask to create a Transitional Layer via the create_barrier_cs function from the leastcostpath package. In the mask Transitional Layer those parts corresponding to the terrain are given a conductance value equal to 1, while everything else (i.e., the parts corresponding to the sea) are given 0 conductance. The mask Transitional Layer is then internally multiplied by the conductance transitional layer representing the cost of movement (according to the user-selected function). This will set to 0 the conductance values of those parts of the study area that do not correspond to the terrain, while keeping unaltered the conductance of those parts that do coincide with the terrain. The irregular.dtm parameter is available in all the implemented functions, movealloc() excluded. Task 17. Can an area where movement is inhibited (barrier) be included in the analysis? Yes. Areas where the movement is inhibited can be fed into the analysis via the barrier parameter; SpatialLineDataFrame or SpatialPolygonDataFrame can be used. The barrier is assigned a conductance value of 0 (i.e., movement is inhibited) by default, but the user can assign any other value via the field parameter. Internally, the barrier creation rests on the create_barrier_cs function from the leastcostpath package. Let’s consider the following example: result1 <- movecost(volc, destin.loc[1,], destin.loc[4,]) result2 <- movecost(volc, destin.loc[3,], destin.loc[6,], move=8) result3 <- movecost(volc, destin.loc[3,], destin.loc[6,], move=8, barrier=result1$LCPs, plot.barrier=TRUE) In the above examples, we use result1 to calculate a LCP that will be used as barrier later on. In result2 , we calculate the LCP between two locations not using any barrier, while in result3 the LCP from result1 is used as barrier. It is apparent that, compared to result2 , the LCP in result3 has to make a long detour to not cross the barrier; the latter is diplayed in blue. If the user does not want the barrier to be displayed, it suffices to set the plot.barrier to FALSE . Note that the move parameter has been set to 8 ; if set to 16 (default), the LCP will be “able” to jump the barrier. Task 18. Can I calculate sub-optimal LCPs between two locations? Yes. The first five sub-optimal LCPs between two locations can be calculated using the moverank() function. See the following example: result <- moverank(volc, destin.loc[3,], destin.loc[6,], funct="t", move=8, lcp.n = 3) In the above plot, the optimal LCP is given 1 (solid black line), while the first two sub-optimal ones are rendered with a different line style. The displayed LCPs are 3 in total, as controlled by the lcp.n parameter that has been set to 3 in this example. The underlying idea is the following: given two locations, we can calculate the least-costly path between them; but, if we disregard that LCP, what path would be the second least costly? And if we in turn disregard those first two, what the third least costly path would be? The same reasoning holds for all the subsequent n-th LCPs. The maximum number of sub-optimal LCPs is currently set at 5. It is worth noting that it may happen that some LCPs will cross one another; this cannot be anticipated and is context dependent. In those cases, the user may want to set the move parameter to 8 (see Task 17 above). As said, that behaviour is context dependent; as a case in point, let’s have a look at another example, which uses in-built data (with a DTM acquired online) and the move parameter set to 16 : result <- moverank(origin=Etna_end_location[3,], destin=Etna_end_location[2,], studyplot=Etna_boundary,funct="wc s", move=16, lcp.n = 3, use.corr = T) In the above example, the wheeled-vehicle critical slope cost function has been used ( funct="wcs" ) and in the resulting plot only the optimal LCP and the first two sub-optimal ones have been rendered. It is worth noting that the first LCP and the third sub-optimal (which runs “next” to it) do not overlap nor cross each other. Task 19. Can the sub-optimal LCPs between two locations be plotted over the leastcost corridor raster? Yes. To achieve that, the user can set the use.corr parameter to TRUE , as we did in the previous task. See the following example: result <- moverank(volc, destin.loc[1,], destin.loc[4,], funct="t", move=8, lcp.n = 3, use.corr=TRUE) In the above plot, the LCPs are plotted on the least-cost corridor between the two locations. It is apparent that the LCPs traverse areas where the cost is the least overall, but only the optimal one (i.e., LCPs n 1) is the one associated to the least cost (see Task 21 below). Task 20. Can barrier(s) be used when calculating sub-optimal LCPs? Yes. To achieve that, the user can feed the barrier into the barrier parameter. See the following example: result1 <- movecost(volc, destin.loc[9,], destin.loc[2,], move=8) result2 <- moverank(volc, destin.loc[1,], destin.loc[4,], move=8, lcp.n = 3, barrier=result1$LCPs, plot.barrier=T RUE) With the above code, we first calculate the LCP between two locations ( result1 ), and then we feed the obtained LCP into moverank() via the barrier parameter; since we want the barrier to be visualised in the rendered plot, we set the plot.barrier parameter to TRUE . As apparent from the above figure, the calcularted LCPs (optimal and sub-optimals) make a detour to avoid the barrier. Needless to say, only three LCPs are rendered since we have set the lcp.n parameter to 3 . Task 21. Can data on the length and cost of each sub-optimal LCP be graphically visualised? Yes. To achieve that, the user can set the add.chart parameter to TRUE . See the following example: result <- moverank(volc, destin.loc[1,], destin.loc[4,], funct="t", time="m", move=8, lcp.n = 3, add.chart=TRUE) The above bubble plot is also rendered. In the chart, the length of the LCPs is plotted against the rank; the size of the bubbles is proportional to the cost of each LCP. It is apparent that the optimal LCP (n 1) is associated to the lowest cost (about 12 minutes walking time). The other sub-optimal LCPs are progressively longer and more costly. It is worth noting that the fourth sub-optimal LCP is slightly shorter than the third, but more costly. All the sub-optimal LCPs are represented in the plot, even though only three overall (the optimal and the first sub-optimal two) are rendered in the plot displaying the study area and the LCPs (see Task 18 about the lcp.n parameter).