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)
Leave a Reply