Welcome to the session on “Spatial Data Analysis in R”. The purpose of this first part is to reproduce with R the raster analysis that you have done yesterday in QGIS. We will repeat almost all steps from yesterday, except the manual drawing of polygon and polyline shapefiles.

1. Preparations and data inspection

Start by loading the packages we need, and by defining our working directory:

library(rgdal)
library(raster)

dir <- "C:/Users/Max/Desktop/gisclass_2023/Day_05_Spatial_Data_Analysis_in_R/Part01_Raster_Analysis/"

You can import shapefiles with the function readOGR() from the package rgdal. Calling the object will show some basic characteristics, including the Coordinate Reference System (crs). You can access the attribute table with the slot “@data”:

cities     <- readOGR(paste0(dir, "data/cities.shp"), verbose=F)
roads      <- readOGR(paste0(dir, "data/roads.shp"), verbose=F)

cities
## class       : SpatialPointsDataFrame 
## features    : 11 
## extent      : -3092740, -1602428, 2730416, 4044825  (xmin, xmax, ymin, ymax)
## crs         : +proj=aea +lat_0=30 +lon_0=95 +lat_1=15 +lat_2=65 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## variables   : 4
## names       : ID,  POP, RASTERVALU, ident 
## min values  :  0,   12,         -9,     1 
## max values  :  0, 3695,        247,     1
roads
## class       : SpatialLinesDataFrame 
## features    : 8 
## extent      : -3017757, -1872973, 2847100, 3751889  (xmin, xmax, ymin, ymax)
## crs         : +proj=aea +lat_0=30 +lon_0=95 +lat_1=15 +lat_2=65 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## variables   : 2
## names       : id, Ã.dentity 
## min values  :  1,         1 
## max values  :  8,         1
cities@data
ID POP RASTERVALU ident
0 630 83 1
0 3683 247 1
0 968 162 1
0 3695 138 1
0 137 203 1
0 12 -9 1
0 99 158 1
0 550 0 1
0 67 0 1
0 1433 0 1
0 780 0 1 
roads@data
id Ã.dentity
0 1 1
1 1 1
2 3 1
3 4 1
4 5 1
5 6 1
6 7 1
7 8 1

You can import a raster with the function raster() from the package raster. By calling the object, you will retrieve some basic information about the resolution, extent, etc.:

dem <- raster(paste0(dir, "data/dem_tiff.tif"))
dem
## class      : RasterLayer 
## dimensions : 1293, 1755, 2269215  (nrow, ncol, ncell)
## resolution : 800, 800  (x, y)
## extent     : -3039518, -1635518, 2828181, 3862581  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=30 +lon_0=95 +lat_1=15 +lat_2=65 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## source     : dem_tiff.tif 
## names      : dem_tiff 
## values     : -32768, 32767  (min, max)

Plot all data on top of each other. Note this would not work if the data had different CRS.

plot(dem)
plot(cities, add=T, pch=20)
plot(roads, add=T)

2. Calculate Slope, Aspect and Proximity

Let’s start with calculating Aspect and Slope from the DEM. We can use the function terrain() for this.

slope  <- terrain(dem, opt="slope",  neighbors=8, unit="degrees")
aspect <- terrain(dem, opt="aspect", neighbors=8, unit="degrees")

par(mfrow=c(2,2), mai = c(0.2, 0.3, 0.5, 0.2))
plot(dem,    main="original DEM", legend=F)
plot(slope,  main="slope",        legend=T)
plot(aspect, main="aspect",       col=rainbow(360), legend=T)
par(mfrow=c(1,1))

Transform the roads to a raster that has the same properties (extent, resolution, crs) as the DEM:

roads_rasterized <- rasterize(roads, dem)
plot(roads_rasterized)

Calculate Euclidean distance to roads:

proximity <- distance(roads_rasterized)
plot(proximity)

3. Reclassification

Now that we have our three basic layers ready (“slope”, “aspect” and “proximity”), we have to reclassify them before we can join all three to calculate a final suitability score. We can use the function reclassify(), which we have to supply with a reclassification matrix:

slope_matrix <- as.matrix(data.frame(from=c(0,2,6), to=c(2,6,20), becomes=c(10,5,0))) 
slope_matrix
##      from to becomes
## [1,]    0  2      10
## [2,]    2  6       5
## [3,]    6 20       0
slope_recl   <- reclassify(slope, slope_matrix, include.lowest=T)
spplot(slope_recl)

Repeat the procedure for aspect:

aspect_matrix <- as.matrix(data.frame(from=c(0,90,270), to=c(90,270,360), becomes=c(0,10,0))) 
aspect_matrix
##      from  to becomes
## [1,]    0  90       0
## [2,]   90 270      10
## [3,]  270 360       0
aspect_recl   <- reclassify(aspect, aspect_matrix, include.lowest=T)
spplot(aspect_recl)

Repeat the procedure for proximity:

proximity_matrix <- as.matrix(data.frame(from=c(0,20000,50000), to=c(20000,50000,500000), becomes=c(10,5,0))) 
proximity_matrix
##       from    to becomes
## [1,]     0 2e+04      10
## [2,] 20000 5e+04       5
## [3,] 50000 5e+05       0
proximity_recl   <- reclassify(proximity, proximity_matrix, include.lowest=T)
spplot(proximity_recl)

4. Calculate suitability score

We can now just add up all reclassified layers to one suitability layer:

suitability <- slope_recl + aspect_recl + proximity_recl
plot(suitability)

Congratulations, you’re done with part 01 and can now move on to part 02!