The GIS operation known as map overlay combines spatial features from one map layer with the attribute (numerical) properties of another.

This R Markdown Notebook explains the R method over, which provides a consistent way to retrieve indices or attributes from a given spatial object (map layer) at the locations of another. It aims to help Percepcion Remota students working at home to understand R functionalities.

library(raster)
## Loading required package: sp
library(sp)

First, clean the memory:

rm(list=ls())

Let’s add a raster object:

(toy <- raster("toy.tif"))
## class      : RasterLayer 
## dimensions : 30, 30, 900  (nrow, ncol, ncell)
## resolution : 0.003679667, 0.002932667  (x, y)
## extent     : -73.93497, -73.82458, 4.946285, 5.034265  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /Users/ivanlizarazo/Documents/ivan/20200223/toy.tif 
## names      : toy 
## values     : 0, 255  (min, max)

Note that this is a single-band toy raster. Its CRS is 4326. But, we would like to have three bands rather than a single one. Let’s create them:

(toy2 <-  as.integer(toy/2))
## class      : RasterLayer 
## dimensions : 30, 30, 900  (nrow, ncol, ncell)
## resolution : 0.003679667, 0.002932667  (x, y)
## extent     : -73.93497, -73.82458, 4.946285, 5.034265  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 0, 127  (min, max)
(toy3 <-  toy2 + 85)
## class      : RasterLayer 
## dimensions : 30, 30, 900  (nrow, ncol, ncell)
## resolution : 0.003679667, 0.002932667  (x, y)
## extent     : -73.93497, -73.82458, 4.946285, 5.034265  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 85, 212  (min, max)
(new_toy <- stack(toy, toy2, toy3))
## class      : RasterStack 
## dimensions : 30, 30, 900, 3  (nrow, ncol, ncell, nlayers)
## resolution : 0.003679667, 0.002932667  (x, y)
## extent     : -73.93497, -73.82458, 4.946285, 5.034265  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      : toy, layer.1, layer.2 
## min values :   0,       0,      85 
## max values : 255,     127,     212

Give a name to each band:

names(new_toy)  <-  c("b1", "b2", "b3")

Plot a color composition:

plotRGB(new_toy, r=1, g=2, b=3)

Let’s remember that this raster is a toy.

Now, let’s add a vector object:

(polys <- shapefile("./polys/polys.shp"))
## class       : SpatialPolygonsDataFrame 
## features    : 10 
## extent      : -73.93497, -73.82511, 4.946305, 5.034265  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       : id, clase 
## min values  :  1,  crop 
## max values  :  9, water

Note that this is a SPDF (i.e a SpatialPolygonsDataFrame) representing land cover classes. Its CRS is 4326.

Let’s visualize the data:

plot(toy)
plot(polys,  add=TRUE)
invisible(text(coordinates(polys), labels=as.character(polys$clase), cex=0.6))

Let’s assume that we would like to obtain a sample of 300 points located within the land cover polygons.

samples <- spsample(polys, 300, type='regular')

What we got? Let’s find it out:

samples
## class       : SpatialPoints 
## features    : 298 
## extent      : -73.93434, -73.82578, 4.946845, 5.033502  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

How many samples?

(numero <- length(samples))
## [1] 298

Note that ‘samples’ is a SpatialPoints object. We wanted 300 points but we may get a different number of points.

The object ‘samples’ can be converted to a SpatialPointsDataFrame:

(df_samples <- as(samples,"SpatialPointsDataFrame"))
## class       : SpatialPointsDataFrame 
## features    : 298 
## extent      : -73.93434, -73.82578, 4.946845, 5.033502  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 0

It is a good idea to add an ID to identify each point. In addition, we probably need an additional attribute, size, for using it when plotting the data.

df_samples@data = data.frame(ID=1:numero, size=1)
df_samples
## class       : SpatialPointsDataFrame 
## features    : 298 
## extent      : -73.93434, -73.82578, 4.946845, 5.033502  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       :  ID, size 
## min values  :   1,    1 
## max values  : 298,    1

An interesting question is where are the samples located?

plot(toy)
plot(polys,  add=TRUE)
plot(df_samples, pch = 1, cex = (df_samples$size)/4, add=TRUE)

Now, it’s time to try the well-known over function. First, let’s use the samples object (i.e. a SpatialPoints object) to overlay the sampling points with the polys object:

##uncomment the following line to see what happens:
##samples$class <- over(samples, polys)$class

The output of the above chunk reads: Error in validObject(.Object) : invalid class “SpatialPointsDataFrame” object: number of rows in data.frame and SpatialPoints don’t match

What happened? It seems that there are several sampling points which are not located within the land cover polygons. This is not what we expected when using the spsample function, right? But, according to the Irishman, “It’s what it is”.

Seriously, a quick look at the error message suggests that we are using an invalid class SpatialPointsDataFrame object. Let’s try again using the df_samples object (i.e. the SpatialPointsDataFrame object).

df_samples$clase <- over(df_samples, polys)$clase

It seems that we did it right this time.

df_samples
## class       : SpatialPointsDataFrame 
## features    : 298 
## extent      : -73.93434, -73.82578, 4.946845, 5.033502  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 3
## names       :  ID, size, clase 
## min values  :   1,    1,  crop 
## max values  : 298,    1, water

Now, let’s extract raster values corresponding to the samples location:

df <- raster::extract(new_toy, df_samples)
head(df)
##       b1  b2  b3
## [1,] 229 114 199
## [2,] 229 114 199
## [3,] 229 114 199
## [4,] 229 114 199
## [5,] 229 114 199
## [6,] 247 123 208

Now, the spectral profiles:

ms <- aggregate(df, list(df_samples$clase), mean)
# instead of the first column, we use row names
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms
# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('green', 'brown', 'yellow', 'cyan', 'blue')
#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,300), xlim = c(1,3), type='n', xlab="Bands", ylab = "Brightness Value")
# add the different classes
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="Spectral Profile from Toy", font.main = 2)
# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

Hope this notebook is helpful.

In case you run into trouble, please read the raster and sp documentation. Look for examples and you will learn a lot.