setwd("YOUR-PATH")

raw_data <- read.table(file = "data/FS_UTM.csv",
                       header = TRUE,
                       sep = ";",
                       stringsAsFactors = FALSE,
                       quote = "")

## check for issues/problems
## str(raw_data)
## apply(X=raw_data[,4:20], MARGIN = 2, FUN=function(x){x[x>1]})
raw_data$settlement[raw_data$settlement > 1] <- 1

library(sp)
sites <- raw_data
coordinates(sites) <- ~xUTM+yUTM
proj4string(sites) <- CRS("+init=epsg:32634")

Please go to http://srtm.csi.cgiar.org/, download the appropriate SRTM scence(s) for our Romanian research area, unzip them and load them into R using the rgdal package.

Plot the sites onto the raster to see whether the srtm scene covers all the points.

library(rgdal)
srtm <- readGDAL("data/srtm_41_03.tif")

## data/srtm_41_03.tif has GDAL driver GTiff 
## and has 6001 rows and 6001 columns

image(srtm)
points(spTransform(sites, CRS("+init=epsg:4326")))

library(raster)
srtm <- raster("data/srtm_41_03.tif")
plot(srtm)
points(spTransform(sites, CRS("+init=epsg:4326")))

## takes a bit of time
srtm <- projectRaster(srtm, res=90, crs=CRS("+init=epsg:32634"))
plot(srtm)
points(sites)

Ask the manual: where is the difference between crop and mask?

srtm <- crop(x = srtm, y = extent(sites)+20000)
plot(srtm)
points(sites)

Since we do not want to reproject and crop in the subsequent analyses we just export the SRTM.

writeRaster(x = srtm,
            filename = "./results/srtm.tif",
            overwrite = TRUE
            )

Using the raster package we can create a lot of different terrain parameters by just one line of code. The result will be a multi-layer raster object, like you know it from multi-/hyperspectral satellite images. You can work on such scenes like you would work on ordinary vector objects in R.

srtm.tp <- terrain(x = srtm,
                   opt = c("slope",
                           "aspect",
                           "TPI",
                           "TRI",
                           "roughness",
                           "flowdir"),
                   unit = "degrees",
                   neighbors = 8)

plot(srtm.tp)

You should be familiar with slope and aspect but do you know the what the other parameter are? A look in the help offers insights (they are based on Wilson & Gallant (2000)):

• "TRI (Terrain Ruggedness Index) is the mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells.
• TPI (Topographic Position Index) is the difference between the value of a cell and the mean value of its 8 surrounding cells.
• Roughness is the difference between the maximum and the minimum value of a cell and its 8 surrounding cells" (from ?raster::terrain)

The help also explains that we can use focal functions in order the adapt the approaches to our needs. A focal function corresponds to a moving window. It uses a matrix of weights for the neighborhood of the focal cells.

tpiw <- function(x, w=5) {
    m <- matrix(1/(w^2-1), nc=w, nr=w) # mean of the surrounding pixel
    m[ceiling(0.5 * length(m))] <- 0 # set the centre cell 0
    f <- focal(x, m, pad = TRUE) # apply moving window
    x - f
}
tpi15 <- tpiw(x = srtm, w=15)
tpi31 <- tpiw(x = srtm, w=31)
#par(mfrow=c(1,3))
#plot(srtm.tp$tpi)
#plot(tpi15)
#plot(tpi31)

For better visualisation we calculate a hillshade and overlay the TPI raster.

srtm.hs <- hillShade(slope = terrain(srtm, opt="slope"),
                     aspect = terrain(srtm, opt="aspect"),
                     angle= 150, direction = 45, normalize = TRUE)

plot(srtm.hs,
     col = gray.colors(n= 255, start=.2, end=.9),
     legend = FALSE)
plot(tpi31,
     col=colorRampPalette(c("red", "white", "blue"))(255),
     alpha = .5, add=TRUE)

writeRaster(x = tpi31,
            filename = "./results/tpi31.tif",
            overwrite = TRUE
            )
list.files(path = "./results")

## [1] "srtm.tif"  "tpi31.tif"

Task: please, make it nice, i.e. visually appealing and informative.

library(rasterVis)
#srtmTheme <- rasterTheme(region=terrain.colors(200))
#levelplot(srtm, par.settings = srtmTheme)
levelplot(brick(tpi31, tpi15,srtm.tp$tpi), par.settings = RdBuTheme)