Skip to content

Raster data

As with vector data, raster data can be handled with various packages. The most common package is the terra package, successor of the deprecated raster package. terra package is also capable process vector data, but it is not as complex as sf package in this way. 

# install.packages("terra")
library(terra)

Reading raster data - SpatRaster

dem <- rast("data/eudem.tif")
dem

plot the raster data 

plot(dem)
plot the histogram of the raster data 
hist(dem)
get more bins

hist(dem, breaks = 100)

Terrain analysis

function terrain can be used to calculate various terrain attributes from the DEM. See ?terrain for more details.

slope <- terrain(dem, "slope")
aspect <- terrain(dem, "aspect")

Working with multiple layers

Class SpatRaster can handle multiple layers (bands) of raster data. These can be accessed by [[ operator as in lists. Layers can be stacked together simply with c() function, but layers have to have the same extent, resolution and crs.

stack the layers together 

dem <- c(dem, slope, aspect)
plot entire object 
plot(dem)

or all histograms of the slope layer 

hist(dem)

Layers can be named, and accessed by name. 

names(dem)
names(dem) <- c("dem", "slope", "aspect")
dem[["dem"]]

Filtering layers can be done by [[ operator with index including negative idenxing.

drop the first layer: 

dem2 <- dem[[-1]]
get the first layer: 
dem3 <- dem2[[1]]

plot only the slope layer from dem object

plot(dem[[2]])
# or
plot(dem[["slope"]])

Summary statistics

How to get mean, min, max, sum, etc. of the raster data?

mean(dem)
max(dem)

This seems not works as expected. That is because theese functions calculates statistics across all layers, eg. maximum value of all layers at each certain pixel position.

To get quick statistics for each layer, you can use summary function.

summary(dem)

or work with the raw values of the raster data via function values(), but take care of the NA values.

mean(values(dem), na.rm = TRUE)

!!! note matrix

Reclassification

logical reclassification

simply use logical operators to create new raster object with values TRUE and FALSE (1 and 0)

reclassified_dem <- dem[[1]]  > 1000
plot(reclassified_dem)
multiple conditions can be combined with logical operators 
reclassified_dem <- dem[[1]] > 500 & dem[[1]] < 550

replacing values

use of subsetting and logical operators can be used to replace values in the raster data. 

reclassified_dem <- dem[[1]]
reclassified_dem[reclassified_dem > 500] <- NA
plot(reclassified_dem)

Tip

functions replace or ifelse works similar way, but can be used also for more complex operations

reclassified_dem <- replace(dem[[1]], dem[[1]] > 1000, 1)
plot(reclassified_dem)

value reclassification

Standard way to reclassify a raster is use of classify function, which reclassify the raster data with reclassification matrix. The matrix has to have 3 columns, first two are the range of values and the third is the value to replace the values in the range. 

reclass_matrix <- matrix(c(0, 500, 1,    
                           500, 1000, 2,  
                           1000, Inf, 3),  
                         ncol = 3, byrow = TRUE)

reclass_matrix

reclassified_dem <- classify(dem[[1]], reclass_matrix)
plot the reclassified raster 
plot(reclassified_dem)

Map algebra

Map algebra is a way to perform operations on raster data. It is similar to vector data, but the operations are performed on the pixel values. 

dem2 <- dem[[1]] * 2
plot(dem2)

sum(dem)

Tip

map algebra is commonmly used for multispectral data, eg. to calculate vegetation indices NDVI, NDWI, etc.

exercise

Find the suitable areas for vineyards, these areas are on southern slopes (aspect between 135 and 225 degrees) with slope between 5 and 15 degrees, and elevation between 200 and 500 meters. 

dem <- rast("/home/ok/git/orthoptera-atlas/processed_data/eudem.tif")
slope <- terrain(dem, "slope")
aspect <- terrain(dem, "aspect")


# Sklon mezi 5-30% (3–16 stupňů)
suitable_aspect <- (aspect >= 135) & (aspect <= 225)
suitable_slope <- (slope > 4) & (slope < 16)
suitable_elevation <- (dem >= 200) & (dem <= 500)

suitable_area <- suitable_aspect & suitable_slope &  suitable_elevation
If you want to keep all layers in the object, you can stack them together. 

suitable_area <- c(suitable_aspect, suitable_slope,suitable_elevation)
plot(suitable_area)

min(suitable_area)
plot(min(suitable_area))
get the area of the suitable areas

sum(suitable_area)
plot(sum(suitable_area))

writing raster data

writeRaster(suitable_area, "outputs/suitable_area.tif", overwrite = TRUE)