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:
- ASTER DEM tile from NASA EarthData
- 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.