# Spectral signature and its standard deviation using R

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.

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

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")
``````
• 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) May 1, 2021 at 8:31
• I will try to do it using the `pivot_longer` function from `tidyr` package. May 1, 2021 at 8:47