Section 15 Using GRASS algorithms
15.1 Concepts
The GRASS software forms an important part of the history of the development of GIS. In the early days of GIS the software was highly specialised. It was assumed at the time that it would be used only by a small number of highly trained technicians and researchers. GRASS was originally designed by the US Army Core of Engineers (USACE). Very large sums of public money were invested in it’s development. The architecture of GRASS was designed for main frame computers running UNIX. The main frame concept no longer exits, so this architecture became an obstacle to using GRASS more widely. However, the main reason behind stagnation in the development of GRASS involved some complex commercial and political considerations.
https://grass.osgeo.org/documentation/applications/
After the code became an open source project elements from GRASS were implemented in other software. The graphical interface to GRASS began to look dated and was difficult to use for less specialised tasks. However GRASS remained popular in some academic applications. The algorithms designed for terrain and hydrological modelling are extremely robust and efficient.
15.1.1 Why use GRASS today?
GRASS algorithms can be called from both QGIS and R without any need to use the native GRASS GUI: There are several areas in which GRASS algorithms have a real edge over other options. The field in which GRASS still excels is traditional hydrological modelling. The algorithms designed for watershed analysis in GRASS are extremely robust and highly customizable for a wide range of uses.
The other area where GRASS can be useful may be rather surprising. GRASS was designed as a raster processing tool-kit. It had a reputation for being rather poor at processing vector data. This was because GRASS used very strict topological rules. Layers produced when digitising in other software with less strict control on the topology could not be easily imported into GRASS. Inconsistencies and errors in the topology led to the layers being rejected by GRASS. This frustrated potential users. However it did lead to the addition of robust algorithms within GRASS for topological checking and error correction. Clean topologies are vital for network analysis. Therefore GRASS is a very useful tool for preparing layers networks for analysis. It also has some very good algorithms for analysing networks.
15.2 Using GRASS from QGIS
Grass algorithms are available from the QGIS toolkit menu.
EXAMPLE TO COME.
15.3 Using GRASS from R
15.3.1 Network analysis in GRASS
In real life applications there is no need to apply bespoke GIS to problems involving driving times and distances. The Google maps API provides direct access to satnav analysis that can be used on PCs, laptops and phones. However there are situations when a bespoke analysis of a road network may still be useful, particularly when working with data from other countries.
A classic network analysis is the travelling salesperson problem. This involves finding the shortest tour that visits all the points that are included only once. The giscourse server holds data from the global road network.
https://sedac.ciesin.columbia.edu/data/set/groads-global-roads-open-access-v1/data-download
Open street map data can also be loaded into R. Let’s try solving the salesperson problem using both as a comparison.
First we will load all the roads in Dorset from the groads-global layer. A query to set up Dorset roads has already been run on the postgis server so this is all that is needed.
library(giscourse)
library(rgdal)
library(dplyr)
library(sp)
library(mapview)
mapviewOptions(basemaps =c("OpenStreetMap","Esri.WorldStreetMap", "Esri.WorldImagery", "OpenStreetMap.BlackAndWhite", "OpenStreetMap.HOT", "Esri.WorldShadedRelief"))
conn<-connect()
rds<-st_read(conn,"dorsetroads")15.3.2 Getting OSM data using the osmdata library
Open street map data has been placed on the server for speed of access. However it can also be pulled down from on-line sources for any area using the osmdata package in R. This is quite tricky to use due to the esoteric structure and nomenclature of OSM data. This will return the OSM roads and tracks within the Dorset roads bounding box.
library(osmdata)
b_box = st_bbox(rds)
dorset_osm <-opq(b_box) %>%
add_osm_feature(key = "highway") %>%
osmdata_sf() %>%
`[[`("osm_lines")
dorset_osm <-dplyr::select(dorset_osm, osm_id, highway)We need to remove some of the minor features.
#unique(dorset_osm$highway)
main<-unique(dorset_osm$highway)[c(1,2,3,4,5,8,11,12,16,18,19,20,21)]
filter(dorset_osm,highway%in%main) ->dorset_osm215.3.3 Linking to GRASS
One of the drawbacks of using both GRASS and SAGA is the fact that data have to be imported into the native format before any analysis can be performed. This can slow down processing time, as the import and export times have to be factored in. However both GRASS and SAGA have extremely efficient algorithms for working on the data once it has been imported. QGIS and R have both adopted the approach of forming “throwaway” temporary layers to hold data imported to SAGA and GRASS format. Importing to GRASS actually involves the formation of a GRASS data base, location and mapset. These are all specifically GRASS concepts that most users do not need to consider in depth. If the linkage in R is made to a named location the data will be saved in GRASS format for use in another R session, so saving reimportation. The easiest way to keep data permanently in GRASS is to use a different location for each CRS. The term “location”" is rather a misnomer in GRASS as the actual geographical location is not important. Consistency in the use of CRS is all that is required.
This line forms, or connects to, a GRASS data base on the server in epsg 4326. The roads layer is used as a template for the location.
link2GI::linkGRASS7(rds, ver_select = TRUE,gisdbase = "/home/rstudio/grass7",location="epsg4326")Now the layers can be written to GRASS.
writeVECT(SDF = as(dorset_osm2, "Spatial"), vname = "osm_roads",v.in.ogr_flags=c("o", "overwrite"))
writeVECT(SDF = as(rds, "Spatial"), vname = "global_roads", v.in.ogr_flags=c("o", "overwrite"))15.3.3.1 Cleaning the topology
The grass7 package runs GRASS commands directly from R. There is extensive documentation of GRASS commands in the manual. The v.clean function is a powerful tool for fixing many topological issues. To run a network analysis of any type it is essential to remove any breaks.
https://grass.osgeo.org/grass70/manuals/v.clean.html
execGRASS(cmd = "v.clean", input = "osm_roads", output = "osm_roads_clean",
tool = "break", flags = "overwrite")
execGRASS(cmd = "v.clean", input = "global_roads", output = "global_roads_clean",
tool = "break", flags = "overwrite")15.3.3.2 Adding some points for the TSP problem
FOr illustration the emap function can be used to digitise some points.
## Run in the console
points<-emap(table="route_points")writeVECT(SDF = as(points, "Spatial"), vname = "points", v.in.ogr_flags=c("o", "overwrite"))15.3.3.3 Finding the solution
To solve the travelling salesperson problem the points are connected to the network of roads. A small threshold of error is added in order to ensure that points that do not fall exactly on top of a line (which has zero width anyway) are linked to it. Then the tsp solving algorithm v.net.salesman is applied.
execGRASS(cmd = "v.net", input = "osm_roads_clean", output = "points1_con",
points = "points", operation = "connect", threshold = 0.1,
flags = c("overwrite", "c"))
execGRASS(cmd = "v.net.salesman", input = "points1_con",
output = "route1", center_cats = paste0("1-", "6"),
flags = c("overwrite","g"))
execGRASS(cmd = "v.net", input = "global_roads_clean", output = "points2_con",
points = "points", operation = "connect", threshold = 0.3,
flags = c("overwrite", "c"))
execGRASS(cmd = "v.net.salesman", input = "points2_con",
output = "route2", center_cats = paste0("1-", "6"),
flags = c("overwrite","g"))library(rgrass7)
osm_route <-rgrass7::readVECT("route1") %>%st_as_sf() %>% st_geometry()
groads_route <- rgrass7::readVECT("route2") %>%st_as_sf() %>% st_geometry()
points<-rgrass7::readVECT("points") %>%st_as_sf() %>% st_geometry()mapview(osm_route,color="red") + mapview(groads_route,color="black") +mapview(points)To improve the analysis from the perspective of actually finding the best route the algorithm should include some consideration of the cost. This can be achieved by assigning values to diffirent types of road.
15.3.4 Watersheds
Watershed analysis takes a digital elevation model as input and finds flows accross the surface. There are many scaling issues to consider when producing output from such an analysis. If a coarse grained DEM is used some of the details of the actual geomorphology may be missed. There are problems with flat areas and sinks formed when a set of pixels fall below the surroundings.
15.3.5 Load a digital evevation model
At the scale of the county of Dorset a model based on SRTM scale of resolution (65m) will produce contextually useful results, as long as the scale of the DEM is taken into account. Using a higher resolution DEM improves the results, but can take an extremely long time to run.
conn<-connect()
dem<-giscourse::mkdem("Dorset",dist=50000,z=9)
dem[dem<0]<-NA
link2GI::linkGRASS7(dem, ver_select = TRUE,gisdbase = "/home/rstudio/grass7",location="epsg27700")If the algorithm is run with a threshold value two useful layers can be generated. One represents the modelled stream flow. This will not coincide exactly with the rivers and streams on the ground due to the introduction of some artefacts and inacuracies at the level of the input DEM. Actual features may also have been modified by human intervention, so streams do not take the route predicted. However the patterns match in general terms. The algorithm can be applied at a much finer scale using LIDAR data in order to analyse site level characteristics for which no other data sources are available. At a regional level there are watershed maps that can be compared with outputs. However the actual output of a watershed model depends on the threshold value. If the threshold is set to a low value the map has a large number of fine scaled basins. At larger thresholds the basins are aggregated. Finding the boundaries of basins can be useful when considering the destination of run off from agricultural fields, among a wide range of other uses.
link2GI::linkGRASS7(dem, ver_select = TRUE,gisdbase = "/home/rstudio/grass7",location="epsg27700")
writeRAST(as(dem, "SpatialGridDataFrame"),"dem", flags = c("overwrite"))
execGRASS(cmd = "r.watershed", elevation = "dem", basin = "basin_hi",stream="stream_hi",threshold =80, flags = c("overwrite", "b"))
execGRASS(cmd = "r.watershed", elevation = "dem", basin = "basin_low",stream="stream_low",threshold =1600, flags = c("overwrite", "b"))
execGRASS(cmd="r.to.vect",input="basin_low",output="basins",type="area", flags = c("overwrite"))
execGRASS(cmd="r.to.vect",input="stream_hi",output="streams",type="line", flags = c("overwrite"))streams<-rgrass7::readVECT("streams")
basins<-rgrass7::readVECT("basins")mapview(basins,alpha.regions = 0) + mapview(streams)dbDisconnect(conn)
rm(list=ls())
lapply(paste('package:',names(sessionInfo()$otherPkgs),sep=""),detach,character.only=TRUE,unload=TRUE)