9.4 GRASS through rgrass7

The U.S. Army - Construction Engineering Research Laboratory (USA-CERL) created the core of the Geographical Resources Analysis Support System (GRASS) (Table 9.1; Neteler and Mitasova 2008) from 1982 to 1995.Academia continued this work since 1997.Similar to SAGA, GRASS focused on raster processing in the beginning while only later, since GRASS 6.0, adding advanced vector functionality (Bivand, Pebesma, and Gómez-Rubio 2013).

We will introduce rgrass7 with one of the most interesting problems in GIScience - the traveling salesman problem.Suppose a traveling salesman would like to visit 24 customers.Additionally, he would like to start and finish his journey at home which makes a total of 25 locations while covering the shortest distance possible.There is a single best solution to this problem; however, to find it is even for modern computers (mostly) impossible (Longley 2015).In our case, the number of possible solutions correspond to (25 - 1)! / 2, i.e., the factorial of 24 divided by 2 (since we do not differentiate between forward or backward direction).Even if one iteration can be done in a nanosecond, this still corresponds to 9837145 years.Luckily, there are clever, almost optimal solutions which run in a tiny fraction of this inconceivable amount of time.GRASS GIS provides one of these solutions (for more details, see v.net.salesman).In our use case, we would like to find the shortest path between the first 25 bicycle stations (instead of customers) on London’s streets (and we simply assume that the first bike station corresponds to the home of our traveling salesman).

  1. data("cycle_hire", package = "spData")
  2. points = cycle_hire[1:25, ]

Aside from the cycle hire points data, we will need the OpenStreetMap data of London.We download it with the help of the osmdata package (see also Section 7.2).We constrain the download of the street network (in OSM language called “highway”) to the bounding box of the cycle hire data, and attach the corresponding data as an sf-object.osmdata_sf() returns a list with several spatial objects (points, lines, polygons, etc.).Here, we will only keep the line objects.OpenStreetMap objects come with a lot of columns, streets features almost 500.In fact, we are only interested in the geometry column.Nevertheless, we are keeping one attribute column; otherwise, we will run into trouble when trying to provide writeVECT() only with a geometry object (see further below and ?writeVECT for more details).Remember that the geometry column is sticky, hence, even though we are just selecting one attribute, the geometry column will be also returned (see Section 2.2.1).

  1. library(osmdata)
  2. b_box = st_bbox(points)
  3. london_streets = opq(b_box) %>%
  4.   add_osm_feature(key = "highway") %>%
  5.   osmdata_sf() %>%
  6.   `[[`("osm_lines")
  7. london_streets = dplyr::select(london_streets, osm_id)

As a convenience to the reader, one can attach london_streets to the global environment using data("london_streets", package = "spDataLarge").

Now that we have the data, we can go on and initiate a GRASS session, i.e., we have to create a GRASS spatial database.The GRASS spatial database system is based on SQLite.Consequently, different users can easily work on the same project, possibly with different read/write permissions.However, one has to set up this spatial database (also from within R), and users used to a GIS GUI popping up by one click might find this process a bit intimidating in the beginning.First of all, the GRASS database requires its own directory, and contains a location (see the GRASS GIS Database help pages at grass.osgeo.org for further information).The location in turn simply contains the geodata for one project.Within one location, several mapsets can exist and typically refer to different users.PERMANENT is a mandatory mapset and is created automatically.It stores the projection, the spatial extent and the default resolution for raster data.In order to share geographic data with all users of a project, the database owner can add spatial data to the PERMANENT mapset.Please refer to Neteler and Mitasova (2008) and the GRASS GIS quick start for more information on the GRASS spatial database system.

You have to set up a location and a mapset if you want to use GRASS from within R.First of all, we need to find out if and where GRASS 7 is installed on the computer.

  1. library(link2GI)
  2. link = findGRASS()

link is a data.frame which contains in its rows the GRASS 7 installations on your computer.Here, we will use a GRASS 7 installation.If you have not installed GRASS 7 on your computer, we recommend that you do so now.Assuming that we have found a working installation on your computer, we use the corresponding path in initGRASS.Additionally, we specify where to store the spatial database (gisDbase), name the location london, and use the PERMANENT mapset.

  1. library(rgrass7)
  2. # find a GRASS 7 installation, and use the first one
  3. ind = grep("7", link$version)[1]
  4. # next line of code only necessary if we want to use GRASS as installed by 
  5. # OSGeo4W. Among others, this adds some paths to PATH, which are also needed
  6. # for running GRASS.
  7. link2GI::paramGRASSw(link[ind, ])
  8. grass_path = 
  9.   ifelse(test = !is.null(link$installation_type) && 
  10.            link$installation_type[ind] == "osgeo4W",
  11.          yes = file.path(link$instDir[ind], "apps/grass", link$version[ind]),
  12.          no = link$instDir)
  13. initGRASS(gisBase = grass_path,
  14.           # home parameter necessary under UNIX-based systems
  15.           home = tempdir(),
  16.           gisDbase = tempdir(), location = "london", 
  17.           mapset = "PERMANENT", override = TRUE)

Subsequently, we define the projection, the extent and the resolution.

  1. execGRASS("g.proj", flags = c("c", "quiet"), 
  2.           proj4 = st_crs(london_streets)$proj4string)
  3. b_box = st_bbox(london_streets) 
  4. execGRASS("g.region", flags = c("quiet"), 
  5.           n = as.character(b_box["ymax"]), s = as.character(b_box["ymin"]), 
  6.           e = as.character(b_box["xmax"]), w = as.character(b_box["xmin"]), 
  7.           res = "1")

Once you are familiar with how to set up the GRASS environment, it becomes tedious to do so over and over again.Luckily, linkGRASS7() of the link2GI packages lets you do it with one line of code.The only thing you need to provide is a spatial object which determines the projection and the extent of the spatial database.First, linkGRASS7() finds all GRASS installations on your computer.Since we have set ver_select to TRUE, we can interactively choose one of the found GRASS-installations.If there is just one installation, the linkGRASS7() automatically chooses this one.Second, linkGRASS7() establishes a connection to GRASS 7.

  1. link2GI::linkGRASS7(london_streets, ver_select = TRUE)

Before we can use GRASS geoalgorithms, we need to add data to GRASS’s spatial database.Luckily, the convenience function writeVECT() does this for us.(Use writeRast() in the case of raster data.)In our case we add the street and cycle hire point data while using only the first attribute column, and name them also london_streets and points.Note that we are converting the sf-objects into objects of class Spatial.In time, rgrass7 will also work with *sf-objects.

  1. writeVECT(SDF = as(london_streets, "Spatial"), vname = "london_streets")
  2. writeVECT(SDF = as(points[, 1], "Spatial"), vname = "points")

To perform our network analysis, we need a topological clean street network.GRASS’s v.clean takes care of the removal of duplicates, small angles and dangles, among others.Here, we break lines at each intersection to ensure that the subsequent routing algorithm can actually turn right or left at an intersection, and save the output in a GRASS object named streets_clean.It is likely that a few of our cycling station points will not lie exactly on a street segment.However, to find the shortest route between them, we need to connect them to the nearest streets segment.v.net’s connect-operator does exactly this.We save its output in streets_points_con.

  1. # clean street network
  2. execGRASS(cmd = "v.clean", input = "london_streets", output = "streets_clean",
  3.           tool = "break", flags = "overwrite")
  4. # connect points with street network
  5. execGRASS(cmd = "v.net", input = "streets_clean", output = "streets_points_con", 
  6.           points = "points", operation = "connect", threshold = 0.001,
  7.           flags = c("overwrite", "c"))

The resulting clean dataset serves as input for the v.net.salesman-algorithm, which finally finds the shortest route between all cycle hire stations.center_cats requires a numeric range as input.This range represents the points for which a shortest route should be calculated.Since we would like to calculate the route for all cycle stations, we set it to 1-25.To access the GRASS help page of the traveling salesman algorithm, run execGRASS("g.manual", entry = "v.net.salesman").

  1. execGRASS(cmd = "v.net.salesman", input = "streets_points_con",
  2.           output = "shortest_route", center_cats = paste0("1-", nrow(points)),
  3.           flags = c("overwrite"))

To visualize our result, we import the output layer into R, convert it into an sf-object keeping only the geometry, and visualize it with the help of the mapview package (Figure 9.3 and Section 8.4).

Shortest route (blue line) between 24 cycle hire stations (blue dots) on the OSM street network of London.
Figure 9.3: Shortest route (blue line) between 24 cycle hire stations (blue dots) on the OSM street network of London.

  1. route = readVECT("shortest_route") %>%
  2.   st_as_sf() %>%
  3.   st_geometry()
  4. mapview::mapview(route, map.types = "OpenStreetMap.BlackAndWhite", lwd = 7) +
  5.   points

Further notes

  • Please note that we have used GRASS’s spatial database (based on SQLite) which allows faster processing.That means we have only exported geographic data at the beginning.Then we created new objects but only imported the final result back into R.To find out which datasets are currently available, run execGRASS("g.list", type = "vector,raster", flags = "p").
  • We could have also accessed an already existing GRASS spatial database from within R.Prior to importing data into R, you might want to perform some (spatial) subsetting.Use v.select and v.extract for vector data.db.select lets you select subsets of the attribute table of a vector layer without returning the corresponding geometry.
  • You can also start R from within a running GRASS session (for more information please refer to Bivand, Pebesma, and Gómez-Rubio 2013 and this wiki).
  • Refer to the excellent GRASS online help or execGRASS("g.manual", flags = "i") for more information on each available GRASS geoalgorithm.
  • If you would like to use GRASS 6 from within R, use the R package spgrass6.