1

I am trying to plot a spectral signature using R similar to the one shown in the image. This plot was obtained using Semiautomatic Classification Plugin-QGIS. Half-transparent regions show the range of values, however I would like to represent the standard deviation.

enter image description here

I extracted the pixel values in each band and its class (land cover). My current data structure is a DataFrame like the following:

red <- c(152, 28, 34, 36, 139, 144, 29, 153)
green <- c(23, 78, 69, 68, 22, 18, 76, 26)
blue <- c(139, 18, 14, 29, 137, 162, 19, 149)
class <- c(1, 2, 2, 2, 1, 1, 2, 1)

df <- data.frame(blue, green, red, class)
df$class <- as.factor(df$class)
names(df) <- c("430nm", "530nm", "620nm", "class")

# my data structure
df

> df
    430   530   620 class
1   139    23   152     1
2    18    78    28     2
3    14    69    34     2
4    29    68    36     2
5   137    22   139     1
6   162    18   144     1
7    19    76    29     2
8   149    26   153     1

I would like to show the standard deviation instead of the range and keep the scale of the X axis and the dashed line based on the wavelength of each band.

1 Answer 1

5

I think that you may want the 95% Confidence Interval and not the Standard Deviation, which is not terribly informative.

The Confidence Interval is the standard error * critical value for the T distribution P( t > cv ) = 0.025 = P( t < -cv )

Your example data was not quite large enough to effectively demonstrate a statistical solution. Sans plotting multiple classes, here is what you are after. The top panel is the Confidence Interval and the bottom is the +/- Standard Deviation.

This creates a transparency color. You can get the RGB values for a named color using col2rgb("blue").

p=75 # pct transparency
clr <- rgb(red=0, green=0, blue=255, max = 255, 
         alpha = (100 - p) * 255 / 100)

Here is some example data. I used a time-series simulation because it yielded a nice serial distribution akin to a spectral signature.

f <- abs(as.vector(arima.sim(list(order = c(1,1,0), ar = 0.7), n = 200)))
  f <- ((f - min(f)) * 255 / (max(f) - min(f)))

Now, plot the standard error and standard deviation. I am using base plot but, you can avoid having to use the polygon function if you created the plot(s) using ggplot2.

par(mfrow=c(2,1)) 

# 95% confidence interval
se <- sd(f)/sqrt(length(f))
cv <- qt(0.95, df=(length(f)-1)) 
L = f -  cv * se
U = f + cv * se 
plot(f, type="l")
  polygon(c(1:length(f), rev(1:length(f))), c(L, rev(U)),  
          col = clr, border = FALSE)
 
# +/- 1 Standard Deviation
L = f -  sd(f)
U = f + sd(f)
plot(f, type="l")
  polygon(c(1:length(f), rev(1:length(f))), c(L, rev(U)),  
          col = clr, border = FALSE)

Here is a multi class example of confidence intervals

d <- rbind(data.frame(class=rep(1,201), w=seq(0,3,0.01)[1:201], 
  dn=abs(as.vector(arima.sim(list(order = c(1,1,0), ar = 0.9), n = 200)))), 
  data.frame(class=rep(2,201), w=seq(0,3,0.01)[1:201], 
  dn=abs(as.vector(arima.sim(list(order = c(1,1,0), ar = 0.6), n = 200)))))

ci <- function(x, p=0.95) qt(p, df=(length(x)-1)) * sd(x)/sqrt(length(x))   
clr <- c(rgb(0, 0, 255, max = 255, alpha = (100 - 80) * 255 / 100),
        rgb(255, 0, 0, max = 255, alpha = (100 - 80) * 255 / 100)) 

# 95% confidence interval
plot(d[d$class == 1,]$w, d[d$class == 1,]$dn, type="n", 
     xlim=range(d$w), ylim=range(d$dn), 
     xlab="digital numbers", ylab="wavelength")

f = d[d$class == 1,]$dn 
L = f - ci(f); U = f + ci(f)
  polygon(c(seq(0,3,0.01)[1:201], rev(seq(0,3,0.01)[1:201])),   
          c(L, rev(U)), col = clr[1], border = FALSE)
    lines(d[d$class == 1,]$w, d[d$class == 1,]$dn)    
f = d[d$class == 2,]$dn 
L = f -  ci(f); U = f + ci(f)
  polygon(c(seq(0,3,0.01)[1:201], rev(seq(0,3,0.01)[1:201])),   
          c(L, rev(U)), col = clr[2], border = FALSE)
    lines(d[d$class == 2,]$w, d[d$class == 2,]$dn)    
abline(v=c(0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2), col="grey")
  legend("topleft", legend=paste0("class", c(1,2)),
         fill = clr, bg="white")
2
  • Thanks for your help. I agree with you. Representing the 95% confidence interval may be a better option than standard deviation, but thank you for providing both options. However my difficulty now is to convert my data frame with 5 columns (430nm 530nm 620nm class) into a dataframe with 3 columns (class, wavelength and DN)
    – sermomon
    May 1, 2021 at 8:31
  • I will try to do it using the pivot_longer function from tidyr package.
    – sermomon
    May 1, 2021 at 8:47

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.