0

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?

Edit - Here is my code (still need help with the for loop to create separate DEMs):

library (rgdal)
library (raster)
#Import the DEM
dem <- raster("Path/to/DEM.tif")
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:
dem.copy <- dem
for (i in elevs) {
  dem.copy[dem.copy>i]=NA
}
  • 1
    you can just use a loop and keep doing raster[raster > 5105] = NA with whatever values you want. values = c(5105, 5225, 5450) for (i in values){raster[raster > i] = NA} – Stefano Potter Nov 20 '18 at 19:51
  • Thanks for the reply. I was able to get the single extraction to work (e.g., raster[raster > 5105] = NA). However, I'm not able to get the multiple raster outputs from the for loop. Do you have any suggestions? TIA! – Lee Walston Nov 20 '18 at 21:58
  • Do you want to save each output condition to a file? – Stefano Potter Nov 20 '18 at 22:09
  • Make a copy of the original raster in the loop, you overwrote with each condition – Stefano Potter Nov 20 '18 at 22:11
  • Hi, yes, I do wish to make each output as a separate file. How do I do that automatically? I wish to have raster5105, raster5225, raster5450, raster5500, and so on... Thanks for your input! :-) – Lee Walston Nov 21 '18 at 2:08
1

You almost had the idea, but since you didn't place your copy inside the loop you already overwrote everything, just put it inside the loop.

library (rgdal)
library (raster)

#Import the DEM   
dem <- raster("Path/to/DEM.tif")
elevs <- c(5175.5, 5176.50, 5177.0, 5177.25, 5178.00)


for (i in elevs){

   #copy raster, important this is inside the loop
   c <- dem

   #mask values
   c[c > i] = NA

   #out_path is whatever you want it to be
   writeRaster(c, out_path)

}
  • I would highly recommend not using R primitives/functions as object names as it can have very undesirable outcomes. – Jeffrey Evans Sep 13 '19 at 14:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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