Fastest Distance Raster Method in R

By: Russell J Gray

Like most people, I have traditionally used the raster package function rasterize() and distance() to generate distance rasters (for example distance to river rasters). However, both of these functions take an incredibly long time, especially when dealing with large, high-resolution data. So here I present a workflow that will take seconds, rather than hours by creating a distance to river raster layer for the entire country of Pakistan as an example:

step 1: download GRASS 7.8 and install on your machine from here:
https://grass.osgeo.org/download/

step 2: run following R code

# open some libraries you will need
library(rgdal)
library(raster)

# set working directory to read in shapefile and country boundary
setwd("")
list.files(pattern=".shp)
water <- readOGR(getwd(), "water_shapefile") # this is the river SpatialLinesDataframe
country <- readOGR(getwd(), "country_shapefile") # this is the country outline

# make an empty raster that matches the extent of your country
r <- raster(x = extent(country), nrows = 2000, ncols = 2000)
r[] <- 1 # set raster value to standard 1 in all cells
# set a latlong projection (this needs to match your shapefiles too)
crs(r) <- "+proj=longlat +datum=WGS84 +no_defs"

# install and open fasterRaster make sure to remove warning errors or it will fail
Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS="true")
remotes::install_github('adamlilith/fasterRaster', dependencies=TRUE)
library(fasterRaster)

# set GRASS directory
grassDir <- 'C:/Program Files/GRASS GIS 7.8'

# make the raster!
distToRiver <- fasterVectToRastDistance(rast=r, # this is the raster
vect=water, # this is the SpatialLinesDataframe
metric="geodesic", # this is set because we are in latlong projection
grassDir=grassDir, # GRASS directory path string
meters = TRUE) # this is also set because we are in latlong projection

# mask the raster to the country shapefile
DistToRiver_m <- mask(DistToRiver, country)
plot(DistToRiver_m, main='Distance to water (m)')
plot(water, col='blue', add=TRUE) 
plot(country, add=TRUE)

And there you have it!