Create raster of aspect relative to the nearest road

Question

I know that slopes which are downhill from the nearest road are an important raster layer for my model (example application: cars carry seeds on their wheels. The seeds fall on the road and downhill from the road, allowing the weeds to spread). It's not enough to have a raster of each cell's distance to the nearest road or a raster of aspect by cardinal direction. I need to know if each cell is downhill or uphill from the nearest road, and if that hill is parallel to verses perpendicular to the nearest road. Ravines which are downhill from, adjacent to, and perpendicular to a nearby road are particularly important. What I'm getting at is I need a raster layer of aspect relative to the nearest road. If my road is perfectly perpendicular to a hill, then to the right of my imaginary car on the road is uphill (arbitrarily assigned zero degrees) and to the left of my car is downhill (180 degrees). The road winds on the hill a little bit, so the entire hill isn't perfectly perpendicular to the road (which is where we get additional degree values). I have multiple roads and multiple hills, so I'll need to be sure to orient each hill to the proper road.

I'm building the raster in R (though it is computationally intense and I may need to switch to Google Earth Engine in the future). For testing, I have the following data for the Mediteranian island of Cyprus:

1. ASTER DEM tile from NASA EarthData
2. road lines shapefile from the Human Data Exchange (which had been imported from OpenStreetMap).

library(raster)
library(rgdal)
library(sp)
library(terra)

dem <- raster::raster('Filepath/ASTGTMV003_N34E032_dem.tif')
roads <- rgdal::readOGR('Filepath/hotosm_cyp_roads_lines_shp/hotosm_cyp_roads_lines.shp') # Note this is a large shapefile and takes a few minutes to read-in

# crop for testing
ext <- extent(32.7, 32.8, 34.9, 35.0)
rds <- crop(roads, ext)
dem <- crop(dem, ext)

# raster::extract() returns the value of the raster at the locations of vector data, so I can get the elevation of the roads at any point, but the extract() function returns a list, not a raster.
el_rd <- extract(x = dem, 
                 y = rds, fun=NULL, na.rm=FALSE, cellnumbers=FALSE, df=FALSE, layer,
                 nl, factors=FALSE, along=FALSE, sp=FALSE)


# I think my workflow will look something like:

# Create new raster of each DEM grid cell to its nearest road

# For each grid cell, compare that elevation value to the elevation of the
# closest point on the closest road (downhill vs uphill from road)

# assign new raster layer based on relative height to and direction of the nearest road

As you can see, I haven't gotten very far. I suspect the raster::extract() and the raster::distance() functions will help. I don't know too much about the terra library, but I've heard it's an improvement to the raster library. I also found this stack overflow page categorizing raster cells by their nearest polygon, but they treat the entire polygon as a single entity with one value, whereas I need each elevation value for every grid cell of the roads. This page calculates vertical distance to the nearest water, but it isn't quite what I need.

I can only work in free open-source platforms. So ArcMap is not an option here.

Answer 1

I am a bit confused here and am wondering if you are overthinking this a bit. Aspect does not indicate uphill/downhill, you are referring to slope. Aspect is the circular direction (0-360 degrees) of the hillslope and slope is the planar angle (0-180).

In wading through your post, best I can tell is that you want a slope raster that is constrained to the roads given a specified distance. You can do this via mask. We will do this via the terra and sf libraries.

Add libraries and set your working directory

library(sf)
library(terra)

setwd("filepath")

Read the elevation (DEM) data and calculate slope in degrees

dem <- rast('ASTGTMV003_N34E032_dem.tif')
  slp <- terrain(dem, "slope")

Read road line data and buffer to 200m

roads <- st_read(hotosm_cyp_roads_lines.shp') 
  b <- st_buffer(roads, 100)

Now, you can create a new slope raster that is constrained to the 200m buffer distance using terra::mask. We nest terra::crop to align the extents. You can write the raster to disk with the terra::writeRaster function.

rd_slp <- mask(crop(slp, ext(b)), vect(b)) 
  plot(rd_slp)