I’ve been having a more sustained play with R and GRASS together. I’ve previously used them, in tandem, for processing satellite images, but haven’t been much further. In this post I’ll look at finding the nearest roads to mountain summits, querying their elevations and presenting summary statistics.

First off, you can run R in GRASS or GRASS in R. I prefer the former as it avoids RAM issues. To do this fire up GRASS, ignore the GUI (for now) and type R at the command line. You can also run RStudio from R by typing rstudio & instead. Amazing! For a full breakdown, see this wiki.

For this work I’ve used the Ordnance Survey Terrain50 elevation model and the OpenRoads dataset. Before you can use the OpenRoads data you’ll need to change the column names, I’ve already solved this problem. Shamefully, I downloaded a Munro dataset some years ago and have no idea where from…

You can import your vector and raster data (roads, summits and elevation model) into GRASS using the GUI, this is probably the easiest method. Following this, you’ll need to filter the OpenRoads data to remove “Not Classified” roads, as these are generally not public and usually dirt tracks (769633 of them in GB). Do this with:

v.extract input=roads@network where="class != 'Not Classified'" output=roads_public

Note, any call using system can be executed at the GRASS console or in terminal, without R. The above has a lot of quotes, so is easiest run in the console or terminal!

Next we can run the following code for our analysis (the v.out.ascii and following cleaning calls would likely be better in readVECT):

library(rgrass7)
library(rgdal)

# Add new columns to vector
system("v.db.addcolumn map=Munros@munro columns='dist INT,to_x INT,to_y INT'")

# Get distance and coordinates of nearest road point
system("v.distance from=Munros@munro from_type=point to=roads_public@munro to_type=line upload=dist,to_x,to_y column=dist,to_x,to_y")

# Read results into R
Munro = execGRASS("v.out.ascii", parameters=list(input="Munros@munro", type="point", columns="NAME,HT,dist,to_x,to_y"), intern=T)

# Clean results
Munro = strsplit(Munro, split="|", fixed=T)
Munro = do.call("rbind.data.frame", Munro)[, -3]
colnames(Munro) = c("Easting", "Northing", "Name", "Height", "Road_dist", "Road_x", "Road_y")
Munro[, -3] = apply(Munro[, -3], 2, as.numeric)
Munro$Easting = round(Munro$Easting)
Munro$Northing = round(Munro$Northing)

# Write data to csv
# Awkward, but easiest way to get point file of road locations!
write.csv(Munro, "~/Munro.csv", row.names=F)

# Make point file of nearest road locations
system("v.in.ascii --overwrite input=/home/melville/Munro.csv output=temp separator=comma x=6 y=7 skip=1")

# Interogate dtm at nearest road points
# Annoyingly, this writes to another temp file
system("r.what --overwrite map=Terrain50@PERMANENT points=temp@munro output=temp.csv separator=comma")

# Read temp file of road heights and extract elevations
x = read.csv("~/temp.csv", header=F)
Munro$Road_height = round(x$V4)

# Calculate some statistics
Munro$Diff = Munro$Height - Munro$Road_height
Munro$Gradient = round(100 * Munro$Diff / Munro$Road_dist)

# Make columns integer to avoid decimal places when exported to a shp file
Munro = Munro[, c(3:5, 8:10)]
for(i in 2:ncol(Munro)){
   class(Munro[, i]) = "integer"
}

# Convert munro summit and nearest road points into lines
x = lapply(1:nrow(Munro), function(i){
   rbind(paste(Munro[i, c(1, 2)], collapse=","),
         paste(Munro[i, c(6, 7)], collapse=","),
         paste("NaN,NaN"))
})
x = do.call("rbind.data.frame", x)
write.table(x,
            "~/temp.txt",
            col.names=F,
            row.names=F,
            quote=F)
system("v.in.lines input=/home/melville/temp.txt output=Munro_lines separator=comma")

# Read lines from GRASS
l = readVECT("Munro_lines")
# Add data frame to lines
l@data = Munro

# Write shp file
writeOGR(l, ".", "Munro", "ESRI Shapefile", overwrite_layer=T)

# Clean up
system("rm ~/Munro*")
system("rm ~/temp*")

Now we’ve a dataset called Munro in our R environment, we can make some summary plots (code below):

WordPress helpfully tiled these plots for me, but they were created with:

png("~/Cloud/Michael/websites/blog/Munro_height.png")
hist(Munro$Height, main="Munro elevation", xlab="Elevation (m)")
dev.off()

png("~/Cloud/Michael/websites/blog/Munro_dist.png")
hist(Munro$Road_dist, main="Munro distance from road", xlab="Distance (m)")
dev.off()

png("~/Cloud/Michael/websites/blog/Munro_road_ht.png")
hist(Munro$Road_height, main="Nearest road elevation", xlab="Elevation (m)")
dev.off()

png("~/Cloud/Michael/websites/blog/Munro_diff.png")
hist(Munro$Diff, main="Difference between Munro and road elevation", xlab="Elevation (m)")
dev.off()

png("~/Cloud/Michael/websites/blog/Munro_grad.png")
hist(Munro$Gradient, main="Average gradient between road and Munro", xlab="Gradient (%)")
dev.off()

png("~/Cloud/Michael/websites/blog/Munro_diff_dist.png")
plot(Diff ~ Road_dist, Munro, main="Elevation difference and distance from road for Munros", xlab="Distance to road (m)", ylab="Elevation difference (m)")
dev.off()

We can also find the “easiest” hills and those furthest from a road:

# Easiest hills
head(Munro[order(Munro$Diff), c(1, 5)])
head(Munro[order(Munro$Gradient), c(1, 6)])

# Furthest from a road
head(Munro[order(Munro$Road_dist, decreasing=T), c(1, 3)])

Finally, I used the Carto plugin for QGIS to upload the dataset to Carto. You can view and interact with it on my Carto map.

View and interact with the live version: https://mikerspencer.carto.com/viz/459140d8-1968-11e7-80ad-0e233c30368f/public_map