How can I create a multi-dimensional spectral plot using R?

I have a raster stack of several layers such as NDVI, Slope, texture, and TCT bands, and would like to create a multi-dimensional spectral signature plot from my ground thruth points representing different vegetation classes.

How can I create a multi-dimensional spectral plot using R based on my raster stack and my ground truth data?

Here is my data:

``````Stack<-stack(NDVI, NDVI_glcm\$glcm_variance, Mean_slope, Sr_TCT)
str(GT_corr)
'data.frame':   38098 obs. of  3 variables:
\$ long  : num  237481 237481 237481 237761 237761 ...
\$ lat   : num  1212239 1212239 1212239 1212240 1212240 ...
\$ M.Info: Factor w/ 23 levels "Bas fond","Bas fond brulee",..: 18 18 18 19 19 19 19 19 20 20 ...
``````

For example, the plot could look like something like this if I would consider only 3 of the layers from my stack:

where each point would represent a ground truth point for a given vegetation class, and a different color would be used for each class.

• Do you have in mind something like n-D plots in ENVI? Nov 27, 2017 at 18:27
– M514
Nov 27, 2017 at 18:34
• ENVI is not a part of ArcGIS. You can use ggplot with facets to compare GT between different layers, assigning color and shape by class. But, without example data is difficult to reproduce an helpful code Nov 27, 2017 at 18:38
• Sorry I thought ENVI was an Harris Geospatial's software implemented in ArcGis. I have now aded a sample code of what my data looks like.
– M514
Nov 27, 2017 at 19:01
• You may want to look at a formal separability metric as well. Here is an R function that calculates 6 separability metrics and plots the class separability as well. rdocumentation.org/packages/spatialEco/versions/0.1-7/topics/… Dec 1, 2017 at 19:32

First, you post data structure... Is useless to make examples.

Here I post a reproducible example with same variables, taking account 3 classes:

``````library(RStoolbox)
library(glcm)
library(raster)

# reproducible example

data(lsat)
data(srtm)

NDVI <- spectralIndices(lsat,red=3,nir=4,indices='ndvi')
NDVI_glcm <- glcm(NDVI)
Mean_slope <- terrain(srtm,unit='tangent')
Sr_Guinea_TCT <- tasseledCap(lsat[[c(1:5,7)]],sat='Landsat5TM')

Guinea_stack <- stack(NDVI,NDVI_glcm\$glcm_variance,Mean_slope,Sr_Guinea_TCT)

GT_corr <- structure(list(x = c(620130, 620130, 621690, 620220, 623880,
627330, 627540, 619800, 619590, 622980, 621480, 623550, 624810,
626670, 627690), y = c(-414450, -417690, -416670, -411360, -410760,
-411000, -411420, -411780, -410880, -418890, -412650, -413250,
-414660, -415380, -415320), class = structure(c(1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L), .Label = c("1",
"2", "3"), class = "factor")), .Names = c("x", "y", "class"), row.names = c(NA,
-15L), class = "data.frame")
coordinates(GT_corr) <- ~x+y
proj4string(GT_corr) <- Guinea_stack@crs
``````

With this information, proceed to extract response from raster stack:

``````# extract response
df <- extract(Guinea_stack,GT_corr,df=T)

# replace object ID with class name/ID
df\$ID <- GT_corr@data\$class
``````

Finally, use ggplot with facets to create a comparison plot, you can use any `geom_*`, I prefer boxplots for this kind of comparison:

``````library(magrittr)
library(tidyr)
library(ggplot2)
library(cowplot) # only for plot style

df %>% gather(Variable,Value,-ID) %>%
ggplot(aes(ID,Value,color=ID))+ geom_boxplot() +
facet_wrap(~Variable,scales = 'free') + xlab('Class')+
scale_color_manual(NULL,values = c('1'='black','2'='red','3'='blue'))
``````

For a n-d plot, use `plot3D` package:

``````library(plot3D)

scatter3D(x=df[,2],y=df[,3],z=df[,4],colvar=NULL,
col = c(rep(1,length.out=5),rep(2,length.out=5),rep(4,length.out=5)),
xlab=names(df)[2],ylab=names(df)[3],zlab=names(df)[4])
legend("topleft",title = 'Class',legend=c("1","2","3"),pch=1,col = c(1,2,4),bg="white")
``````

Case example:

Using provided data:

``````library(plot3D)

df <- df[,-1] #delete X column (Feature ID from exported data.frame... Not necessary)
``````

First, identify which classes are stored en ID field:

``````classes <- unique(df\$ID)
``````

Assing one color per class... In this case I used rainbow color palette, but you can use a customized one:

``````col_class <- rainbow(n=length(classes)) # Also, you can set color by name (i.e. 'red') or number (i.e. 10)
df\$color <- col_class[match(df\$ID,classes)]
``````

Plot `scatter3D` and add legend:

``````scatter3D(x=df[,2],y=df[,3],z=df[,4],colvar=NULL, col = df\$color, pch=20,
xlab=names(df)[2],ylab=names(df)[3],zlab=names(df)[4])
legend("topleft",title = 'Class',legend=classes,pch=20,
cex=0.8,y.intersp=1,col = col_class,bg="white")
``````

The link provided above shows several ways to use this package. One of interest is fancy plot:

``````scatter3D_fancy <- function(x, y, z,..., colvar = z)
{
panelfirst <- function(pmat) {
XY <- trans3D(x, y, z = rep(min(z), length(z)), pmat = pmat)
scatter2D(XY\$x, XY\$y, colvar = colvar, pch = ".",
cex = 2, add = TRUE, colkey = FALSE)

XY <- trans3D(x = rep(min(x), length(x)), y, z, pmat = pmat)
scatter2D(XY\$x, XY\$y, colvar = colvar, pch = ".",
cex = 2, add = TRUE, colkey = FALSE)
}
scatter3D(x, y, z, ..., colvar = colvar, panel.first=panelfirst,
colkey = list(length = 0.5, width = 0.5, cex.clab = 0.75))
}

scatter3D_fancy(x=df[,2],y=df[,3],z=df[,4],colvar=NULL, col = df\$color,
xlab=names(df)[2],ylab=names(df)[3],zlab=names(df)[4])
legend("topleft",title = 'Class',legend=classes,pch=1,
cex=0.5,y.intersp=1,col = col_class,bg="white")
``````

• Thank you for this answer! Would it be possible to have one single plot where each band is an axis? Where one could see in an n-dimensional space (here n= 6) the distribution of my ground truth points along the combination of their extract value for each dimension (i.e. layer) of the raster stack? Don't know if I'm clear here..
– M514
Nov 27, 2017 at 20:42
• I have added an example of what I'm trying to do (here in a 3D space). I think I found a package ("plot3D") in which I could do so!
– M514
Nov 27, 2017 at 20:51
• I have tried the code and it works perfectly - thank you very much. How can I change the number of classes here? Let's say from n=3 to n=23?
– M514
Nov 28, 2017 at 20:23
• You need to have a color palette associated with your classes. You can create a new column in output data.frame from `extract()` defining color per observation. Is quite simple using `unique()`, `match()` and desired colors. Save a sample of `extract()` output data.frame, upload it to wetransfer or other service and I'll show you how to do it Nov 28, 2017 at 21:28
• This is great, exactly what I was looking for, Thanks again!!
– M514
Nov 29, 2017 at 16:03