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).

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.

  • 1
    In the future, please do not put the onerous on us to search for data. You just aimed us at the source for the DEM with no other information. If you are going to provide data, package it up and provide a download location. We should not have to do extra work in trying to help you! Apr 8, 2022 at 21:28
  • Google Earth Engine is not an open-source platform, so that is not an option either.
    – Spacedman
    Apr 10, 2022 at 10:16
  • It sounds like you want to do a conventional hydrological flow model, but where the "rain" falls only on the roads. I don't see anything in your outline plan that might cope with a grid cell where the nearest road is over the other side of a hill, such that nothing from that road could "flow" to the grid cell.
    – Spacedman
    Apr 10, 2022 at 10:22
  • @Spacedman Thank you for the recommendation of treating this like a hydrological process. Yes, that is what I'm going for, and I hadn't thought of it that way. It is a weakness that if my grid cell is close to a road horizontally, but up and then down a hill, I haven't accounted for it. Based on your suggestion, I think I need three raster layers: 1) distance to the nearest road (I've already made this) 2) vertical distance to the nearest road - could be negative it the nearest road is lower than my cell. 3) distance to the nearest road which has a higher elevation than my grid cell. Apr 10, 2022 at 17:15
  • I could use Google Earth Engine. I have a GEE account, though I'm not as familiar with that platform. I've started searching for hydrology packages in R, which will hopefully lead me in the right direction. I've found a resource to help me: hess.copernicus.org/articles/23/2939/2019/hess-23-2939-2019.pdf Apr 10, 2022 at 17:17

1 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



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)) 
  • Thank you for the feedback. No, I am not looking to mask or crop based on distance to road. I am trying to create a layer based on slope orientation to road. A description of the raster layer I want in more detail: Let's say I'm driving on a road that is perpendicular to a hill. To the right of my car is uphill (given an arbitrary value of 0 degrees), to the left of my car is downhill (180 degrees). My car is dropping seeds from its wheels. The seeds can only travel downhill. So I need to know which slopes are downhill from my road. I'm editing my post to me more clear. Apr 9, 2022 at 1:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.