2
 r<-raster("DEM.tif")

How to do the classification based on every ten meters for example and have different colors too?

6

You can use the reclassify function in the raster package to reclassify the DEM. The general idea is to generate a reclass matrix which provides the instructions on how to reclassify the continuous DEM elevation values.


enter image description here

require(rgdal)
require(raster)

# Read DEM and convert to raster layer object
dem = raster("C:/temp/dem.tif")

# Generate a reclass matrix (begin, end, value ...)
m = c(290, 300, 1, 300, 310, 2, 310, 320, 3, 320, 330, 4, 330, 340, 5)
rclmat = matrix(m, ncol=3, byrow=TRUE)

# Reclassify the raster layer
reclassified_dem = reclassify(dem, rclmat)

# plot the results
par(mfrow=c(1,2))
plot(dem)
plot(reclassified_dem)

# Write to disk
writeRaster(reclassified_dem, "C:/temp/out_dem.tif", drivername = "GTiff", type = "Float32")
| improve this answer | |
  • very nice,since i'm a beginner explain how you managed to plot 2 things without adding add=T in the second? I guess something with the par happens? Also tell me how you zoom in at a certain point of the map in both the split images – gsa Apr 28 '15 at 13:59
  • par is used to set graphic properties. In this case, I specified that I want to plot two figures in a 1 row X 2 column arrangement. The map is not zoomed in--it is the entire extent of my DEM. – Aaron Apr 28 '15 at 14:37
  • the "write to disk line" what exactly does?exports it? or maybe takes it from memory to disk?if yes what does it mean practically? – gsa Apr 28 '15 at 14:47
  • Note that I have altered the script to make it easier to understand. The # Write to disk line takes a RasterLayer class (i.e. reclassified_dem) which is available only in memory and writes it to your hard drive as a geotiff. You can then use the geotiff in your choice of GISs. – Aaron Apr 28 '15 at 14:56
  • Therefore exports it correct?Why didn't we set location directory for the saved file? – gsa Apr 28 '15 at 15:00
3

This is an alternative using the function cut to assign elevation values into classes of elevation, making possible to discretize colors in the map (color per class of elevation).

#Generate reproducible example
library(raster)
f = system.file("external/test.grd", package="raster") #path to raster file
DEM = raster(f) #import raster file
DEM = rasterToPoints(DEM) #change raster to matrix object with coordinates
DEM = data.frame(DEM) #change matrix to data.frame object

#assign the DEM elevation values to class intervals
DEM$test = cut(DEM$test,breaks=c(100,200,300,400,500,600,700,800,900,1500)) #here you'll want to create classes with intervals of 10 meters.

#Plot raster
library(ggplot2)

ggplot() +
geom_tile(data=DEM,aes(x=x,y=y,fill=test)) + 
  scale_fill_brewer("Elevation (m)", type = "seq", palette = "Greys") +
  xlab("Meters") + ylab("Meters") +
  coord_equal() + 
  theme_bw() + 
  theme(panel.border=element_rect(color="black"),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        axis.title.y=element_text(size = 14, vjust=+0.7),
        axis.title.x=element_text(size = 14, vjust=-0.2),
        axis.text.y=element_text(size = 11),
        axis.text.x=element_text(size = 11),
        legend.title=element_text(size = 13),
        legend.text=element_text(size = 13))

enter image description here

| improve this answer | |
  • +1 Nice approach and it is always nice to see ggplot figures! – Aaron Apr 28 '15 at 14:32

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.