Exercises With Raster Data: Solutions (Part 1 and 2)

Below are the solutions to these exercises on Raster Data Analysis (Part 1 and 2).

####################
#                  #
#    Exercise 1    #
#                  #
####################

library(raster)

## Create a folder named data-raw inside the working directory to place downloaded data     
if(!dir.exists("./data-raw")) dir.create("./data-raw")

## If you run into download problems try changing: method = "wget"        
download.file("https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/srtm_pnpg.zip", 
              "./data-raw/srtm_pnpg.zip", method = "auto")

## Uncompress the zip file        
unzip("./data-raw/srtm_pnpg.zip", exdir = "./data-raw")

# Load raster file into R
rst <- raster("./data-raw/srtm_pnpg.tif")

# Number of rows, columns and layers        
dim(rst)

## [1] 579 555   1

# Number of rows
nrow(rst)

## [1] 579

# Number of columns
ncol(rst)

## [1] 555

# Number of layers
nlayers(rst)

## [1] 1

# Number of cells
ncell(rst)

## [1] 321345

####################
#                  #
#    Exercise 2    #
#                  #
####################

# Spatial resolution (pixel size in x and y dimensions)
res(rst)

## [1] 80 80

# Coordinate reference system: Datum/Ellipsoid: WGS 1984, Projection: UTM 29N, Units: meters
crs(rst)

## CRS arguments:
##  +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0

####################
#                  #
#    Exercise 3    #
#                  #
####################

# The extent object
extent(rst)

## class       : Extent 
## xmin        : 549619.7 
## xmax        : 594019.7 
## ymin        : 4613377 
## ymax        : 4659697

# Lenght in meters
xmax(rst) - xmin(rst)

## [1] 44400

# Height in meters
ymax(rst) - ymin(rst)

## [1] 46320

####################
#                  #
#    Exercise 4    #
#                  #
####################

# Mean
cellStats(rst, mean)

## [1] 747.2759

# Standard-deviation
cellStats(rst, sd)

## [1] 311.8615

####################
#                  #
#    Exercise 5    #
#                  #
####################

# Raster quantiles
cellStats(rst, stat = function(x, ...) quantile(x, probs = c(0.01, 0.25, 0.5, 0.75, 0.99), ...))

##   1%  25%  50%  75%  99% 
##   72  527  774  983 1340

####################
#                  #
#    Exercise 6    #
#                  #
####################

# Get values for the raster (memory may be an issue)
x <- values(rst)

# Make the plot (notice higher deviations for values above the average)
qqnorm(x)
qqline(x)

####################
#                  #
#    Exercise 7    #
#                  #
####################

set.seed(12345)

# Generate random points with uniform distribution bounded by the raster extent
xyRandPoints <- data.frame(x = runif(100, xmin(rst), xmax(rst)), y = runif(100, ymin(rst), ymax(rst)))

# Convert the initial data frame into a SpatialPoints objects for clarity
xyRandPoints <- SpatialPoints(xyRandPoints, proj4string = crs(rst))

# Extract values
xyVals <- extract(rst, xyRandPoints)

print(xyVals)

##   [1] 1217  820  934  825  317  302  669 1101 1028  899  661  206  689  184
##  [15] 1114  975  874  782  319  906  538  623 1079  560  578 1116  999 1229
##  [29]  456 1039  710  397 1104  531 1183  938  942  971  610  627 1190  721
##  [43] 1057  695  320  924   45  135  143 1139  628  975 1184  332  888  977
##  [57]  548  269  988 1012  540  884  612  781 1185  348  851  488 1159 1037
##  [71] 1021 1029 1112 1102  113  874  883 1157 1118  728 1038  760  484  379
##  [85]  215 1112  788 1150  817  719 1009  788  186 1198 1000  619 1109  766
##  [99]  433  566

####################
#                  #
#    Exercise 8    #
#                  #
####################

rstFt <- rst / 0.3048

rstStack <- stack(rst, rstFt)

print(rstStack)

## class       : RasterStack 
## dimensions  : 579, 555, 321345, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 80, 80  (x, y)
## extent      : 549619.7, 594019.7, 4613377, 4659697  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : srtm_pnpg.1, srtm_pnpg.2 
## min values  :     9.00000,    29.52756 
## max values  :    1520.000,    4986.877

####################
#                  #
#    Exercise 9    #
#                  #
####################

# Bounding coordinates
xmin <- 554615
xmax <- 589015 
ymin <- 4618355 
ymax <- 4654705

# Create the extent object by defining the bounding coordinates
newExtent <- extent(xmin, xmax, ymin, ymax)

# Crop
cropRst <- crop(rst, newExtent)

print(cropRst)

## class       : RasterLayer 
## dimensions  : 455, 430, 195650  (nrow, ncol, ncell)
## resolution  : 80, 80  (x, y)
## extent      : 554579.7, 588979.7, 4618337, 4654737  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : srtm_pnpg 
## values      : 47, 1520  (min, max)

####################
#                  #
#    Exercise 10   #
#                  #
####################

# Target CRS
targetCRS <- CRS("+init=epsg:3035")

# Reproject data
rst_ETRS89_LAEA <- projectRaster(rst, crs = targetCRS, res = 100, method = 'bilinear')

print(rst_ETRS89_LAEA)

## class       : RasterLayer 
## dimensions  : 572, 547, 312884  (nrow, ncol, ncell)
## resolution  : 100, 100  (x, y)
## extent      : 2799171, 2853871, 2238356, 2295556  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : srtm_pnpg 
## values      : 10.25648, 1515.846  (min, max)

Related exercise sets:

Spatial Data Analysis: Introduction to Raster Processing (Part 1) Spatial Data Analysis: Introduction to Raster Processing (Part 2) Exercises With Raster Data (Part 1 and 2) Big Data Manipulation in R Exercises Data Manipulation with data.table (part -2)