4

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:

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 '17 at 18:27
  • Unfortunately I don't have access to ArcGIs.. – M514 Nov 27 '17 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 '17 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 '17 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/… – Jeffrey Evans Dec 1 '17 at 19:32
4

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

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

enter image description here


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

enter image description here

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


 Case example:

Using provided data:

library(plot3D)

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

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

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