Qfyilyi

I'd like to use the Raster package in R to process a single DEM. I'd like to use a pre-determined list of elevation values to extract multiple raster datasets from this DEM. Not all values are of equal interval. For example, my DEM ranges between 5,000 and 6,000 ft. I'd like to use a pre-determined list of 10 values to clip this raster (similar to "extract by attribute" tool in Spatial Analyst). The values do not have equal intervals, for example: 5105 5225 5450 5500 . . . and so on...

At each step, I'd like to extract all values LESS than the specific value (e.g., where Value < 5,105).

If I were to do this manually, I'd set up a batch run to repeatedly use the "extract by attribute" tool in Spatial Analyst. I don't want to do that. I will have a lot of input DEMs to do in the future and I'd like to develop a script to chug through these quickly.

Anyone have any ideas?

Here is my code (still need help with the for loop to create separate DEMs). The code below does not produce separate output rasters at each elevation slice. I need to create separate output rasters (the number of output rasters equal the number of values in the elevs list). TIA!

library (rgdal) 

library (raster) 

#Import the DEM 

dem <- raster("Path/to/DEM.tif") 



#List of Elevations

elevs = c(5175.5, 5176.50, 5177.0, 5177.25, 5178.00) 



#Extract DEM at at elevations less than elevs list 

#This can be done manually as follows: 

dem.5175.5 <- dem    

dem.5175.5[dem.5175.5>5175.5]=NA 



#Trying to do this iteratively through the list of elevs

#Need help here...



dem.copy <- dem 

for (i in elevs) {   

   dem.copy[dem.copy>i]=NA 

} 

I don't have a DEM file to test with, but I would write a helper function and use lapply to produce a list of dems split by your elevs vector. Here is a quick example:

DemSplit <- function(dem, elev){

  dem[dem>elev] <- NA

  dem

}



lapply(elev, function(x) DemSplit(dem, x))