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)
'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:

enter image description here

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?
    – aldo_tapia
    Nov 27, 2017 at 18:27
  • Unfortunately I don't have access to ArcGIs..
    – 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
    – aldo_tapia
    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
  • 1
    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

1 Answer 1


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

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


# reproducible example


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

enter image description here enter image description here

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(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'))

enter image description here

For a n-d plot, use plot3D package:


          col = c(rep(1,length.out=5),rep(2,length.out=5),rep(4,length.out=5)),
legend("topleft",title = 'Class',legend=c("1","2","3"),pch=1,col = c(1,2,4),bg="white")

enter image description here

For more information about this package and set-up visit this blog.

 Case example:

Using provided data:


df <- read.csv('path/to/sample_data_df.csv')

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,
legend("topleft",title = 'Class',legend=classes,pch=20,
       cex=0.8,y.intersp=1,col = col_class,bg="white")

enter image description here

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,
legend("topleft",title = 'Class',legend=classes,pch=1,
       cex=0.5,y.intersp=1,col = col_class,bg="white")

enter image description here

  • 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
    – aldo_tapia
    Nov 28, 2017 at 21:28
  • 1
    This is great, exactly what I was looking for, Thanks again!!
    – M514
    Nov 29, 2017 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.