This is an excercise for Advanced Remote Perceptin (PRA) In this excercice I execute a spatial data analysis with data and code of João Gonçalves.

In the first time we should install library that we need (sp, raster, rgdal)

Install dependeces:

#if(!("rgdal" %in% installed.packages()[,1]))
  #install.packages(c("rgdal"), dependencies = TRUE)

#if(!("sp" %in% installed.packages()[,1]))
  #install.packages(c("sp"), dependencies = TRUE)

#if(!("raster" %in% installed.packages()[,1]))
  #install.packages(c("raster"), dependencies = TRUE)

Next Download example data, data is in SRTM - version 4.1

Download Data:

## 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")

Now we’ll read raster data, for this we’ll use raster function (previously we installed library raster) and then we use print command for see information of the raster:

Call Raster:

library(raster)
## Loading required package: sp
# In this example the function uses a string with the data location
rst <- raster("./data-raw/srtm_pnpg.tif")

# Print raster attributes
print(rst)
## class       : RasterLayer 
## dimensions  : 579, 555, 321345  (nrow, ncol, ncell)
## 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 
## data source : D:\A2ML\Documents\data-raw\srtm_pnpg.tif 
## names       : srtm_pnpg 
## values      : 9, 1520  (min, max)

From the above, we can see some important information about our raster dataset. In this exapmple we “call” the raster in our disk with raster function and then we put this in “rst” variable, note that we can call variab1le with any name.

To check if the raster dataset is stored in RAM we use inMemory function:

inMemory(rst)
## [1] FALSE

The package also provides several functions to access each raster attribute individually.

rst <- raster("./data-raw/srtm_pnpg.tif")
## Raster layer name(s) / more useful for multi-layer rasters
## By default coincides with the file name without extension
names(rst)
## [1] "srtm_pnpg"
## Number of rows, columns and layers
dim(rst)
## [1] 579 555   1
## Nr of rows
nrow(rst)
## [1] 579
# Nr of columns
ncol(rst)
## [1] 555
## Total number of grid cells
ncell(rst)
## [1] 321345
## Spatial resolution in x and y dimensions
res(rst)
## [1] 80 80
## Data type - see ?dataType for more details
dataType(rst)
## [1] "INT2S"
## Extent (returns a S4 object of class "Extent")
extent(rst)
## class       : Extent 
## xmin        : 549619.7 
## xmax        : 594019.7 
## ymin        : 4613377 
## ymax        : 4659697

Info on extent coordinates can be retrieved individually:

xmin(rst)
## [1] 549619.7
xmax(rst)
## [1] 594019.7
ymin(rst)
## [1] 4613377
ymax(rst)
## [1] 4659697
#See Coordenate of reference System:
crs(rst)
## CRS arguments:
##  +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0

Above we can see statics calulated, include coordenate system and x,y max y min.

In this part the author explore visualization functions

To visualize the data we can simply use the function plot

rst <- raster("./data-raw/srtm_pnpg.tif")
plot(rst, main="Elevation (meters) for Peneda-Geres\n National Park, Portugal", 
     xlab="X-coordinates", ylab="Y-coordinates")

Note that plot show a image of the raster, main show a tittle, xlab and ylab show a tittle in the axis, furthermore note that the legend is in “”

The author too shows that We can also use a histogram to visualize the distribution of elevation values in a raster

Generate histogram from a sample of pixels (by default 100K are randomly used)

  hist(rst, col="blue", main = "Histogram of elevation", prob = TRUE, 
     xlab = "Elevation (meters a.s.l.)")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 31% of the raster cells were used. 100000 values used.
# Generate the density plot object and then overlap it
  ds <- density(rst, plot = FALSE)
  lines(ds, col = "red", lwd = 2)

For calculate statics variables we can use command summary

summary(rst)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (31.12% of all cells)
##         srtm_pnpg
## Min.            9
## 1st Qu.       527
## Median        774
## 3rd Qu.       984
## Max.         1520
## NA's            0

We can use individual commad to see this statics

Min

minValue(rst)
## [1] 9

Max

maxValue(rst)
## [1] 1520

The package also provides a more general interface to calculate cell statistics using the cellStats function (no sample employed).

Mean

cellStats(rst, mean)
## [1] 747.2759

Standard-deviation

cellStats(rst, sd)
## [1] 311.8615

Median

cellStats(rst, median)
## [1] 774

Median-absolute deviation (MAD)

cellStats(rst, mad)
## [1] 336.5502

Quantiles

5%, 25%, 50%, 75% and 95%

cellStats(rst, function(x,...) quantile(x, probs=c(0.05, 0.25, 0.5, 0.75, 0.95),...))
##   5%  25%  50%  75%  95% 
##  186  527  774  983 1224

Extracting Values

The raster package allows several possibilities to extract data at specific points, lines or polygons. The extract function used for this purpose allows a two-column matrix or data.frame (with x, y coordinates) or spatial objects from the sp package such as: SpatialPoints, SpatialPolygons, SpatialLines or Extent as input.

For the first example, we will start by extracting raster values using points as input:

One specific point location (with coordinates in the same CRS as the input raster)

xy <- data.frame(x = 570738, y = 4627306) # Add Coordenate to Point
xy <- SpatialPoints(xy, proj4string = crs(rst)) # Add Projection, (Same of Raster)
extract(rst, xy)
##     
## 611

Extract raster values for 20 randomly located points

xy <- data.frame(x = runif(20, xmin(rst), xmax(rst)), y = runif(20, ymin(rst), ymax(rst)))
xy <- SpatialPoints(xy, proj4string = crs(rst))
extract(rst, xy)
##  [1]  820 1141  211 1003  921 1031 1128  761 1155  742  788  829  906 1067
## [15]  395  512  234  896  847  742

Typically, we are also interested in extracting raster values for specific regions-of-interest (ROI). In this example, we will use a polygon (a broad-leaf forest area) to assess the distribution of elevation values within it.

Download the vector data with the woodland patch ROI

If you run into download problems try changing: method = “wget”

#download.file("https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/WOODS_PNPG.zip", "./data-raw/WOODS_PNPG.zip", method = "auto")

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

## Convert the data into SpatialPolygons (discards the attached attribute but keeps geometry)
#woods <- as(readOGR(dsn = "./data-raw", layer = "woods_pnpg"), "SpatialPolygons")
woods <- shapefile("D:\\A2ML\\Documents\\data-raw\\woods_pnpg.shp")
## Warning in rgdal::readOGR(dirname(x), fn, stringsAsFactors =
## stringsAsFactors, : Z-dimension discarded

Let’s check out the polygon data with a simple plot:

#woods
rst
## class       : RasterLayer 
## dimensions  : 579, 555, 321345  (nrow, ncol, ncell)
## 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 
## data source : D:\A2ML\Documents\data-raw\srtm_pnpg.tif 
## names       : srtm_pnpg 
## values      : 9, 1520  (min, max)

Plot elevation raster

plot(rst, main="Elevation (meters) for Peneda-Geres\n National Park, Portugal", 
     xlab="X-coordinates", ylab="Y-coordinates") # Plot the raster with comments
plot(woods, add = TRUE) ## Add the ROI

Now, let’s extract the raster values from the polygon ROI and calculate some statistics:

elev <- extract(rst, woods)[[1]] ## Subset the first (and only) geometry element

fivenum(elev)# Tukey's five number summary:  minimum, lower-hinge, median, upper-hinge, and, maximum
## [1]  556  677  745  822 1086

When using extract with a SpatialPolygons* object, by default, we get a list containing a set of raster values for each individual polygon. Now, using the extracted values, we can investigate the distribution of elevation values for the target patch.

hist(elev, main = "Histogram of ROI elevation", xlab = "Elevation (meters a.s.l.)", prob = TRUE)
abline(v = mean(elev), lwd = 2) ## Mean line