Specim IQ test of grey standards

• Agustin.Lobo@geo3bcn.csic.es
• 20221115

This script extracts reflectance spectra of 2 grey standards from a VISNIR hyperspectral image acquired with a Specim IQ system in HSILab premises on 20221109:

• Specim IQ white reference card.
• Labsphere 50% reflectance standard
• Cubert grey reflectance standard

defaultW <- getOption("warn")
options(warn=-1)
#Paths
RSpectDir       <- "/home/alobo/owncloudRSpect/RSpect"
datadir         <- "/home/alobo/WORK/HSI_Lab/Tests_Hypercams/SpecimIQ/20221109/034"
imadir          <- file.path(datadir,"results")
vecdir          <- file.path(datadir,"QGIS/Polygons")
#libraries
library(reshape2, quietly = TRUE)
library(plyr, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(plotly, quietly = TRUE)
library(terra, quietly = TRUE)
library(tidyterra, quietly = TRUE)
library(gridExtra, quietly = TRUE)
library(directlabels, quietly = TRUE)

#Required information
load(file.path(RSpectDir,"IQbands.rda")) #Bands Specim IQ system

1. Read images and polygons

1.1 RGB image

A conventional RGB image in png format is automatically recorded by a dedicated CCD in Specim IQ. By some reason, this image is rotated vs. the hyperspectral one, so the code rotates it.

options(warn=-1)
r <- rast(file.path(imadir,"RGBBACKGROUND_034.png"))

Warning: [rast] unknown extent

r

class       : SpatRaster 
dimensions  : 645, 645, 4  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 645, 0, 645  (xmin, xmax, ymin, ymax)
coord. ref. :  
source      : RGBBACKGROUND_034.png 
names       : RGBBACK~D_034_1, RGBBACK~D_034_2, RGBBACK~D_034_3, RGBBACK~D_034_4 

#plotRGB(r, stretch="lin")
plotRGB(flip(t(r)))

1.2 Reflectance image

The reflectance image was calculated by the Specim IQ system itself using the average value of the White reference

class       : SpatRaster 
dimensions  : 512, 512, 204  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 512, 0, 512  (xmin, xmax, ymin, ymax)
coord. ref. : Arbitrary 
source      : REFLECTANCE_034.dat 
names       :    397.32,    400.20,    403.09,    405.97,    408.85,    411.74, ... 
min values  : 0.1458333, 0.1694915, 0.1578947, 0.1359223, 0.1205674, 0.1185567, ... 
max values  : 1.1458334, 1.1355932, 1.1052631, 1.1067961, 1.0845070, 1.0773196, ... 

We display both RGB and CIR color composites with the following spectral bands:

RGB

	Band	Wavelength
70	Band070	598.60
53	Band053	548.55
19	Band019	449.35

CIR

	Band	Wavelength
203	Band203	1000.49
70	Band070	598.60
53	Band053	548.55

  crs(s) <- ""
  #suppressMessages(
  ggRad <- ggplot() + 
  geom_spatraster_rgb(data = stretch(subset(s, subset=defbands), minq=0.02, maxq=0.85)) +
    coord_fixed(ratio = 1) +
    theme_void() + theme(plot.title = element_text(hjust = 0.5)) +
    ggtitle(paste0("\nBands ", paste(defbands,collapse=", ")))
  ggRad2 <- ggplot() + 
    geom_spatraster_rgb(data = stretch(subset(s, subset=defbands2), minq=0.02, maxq=0.85)) +
    coord_fixed(ratio = 1) +
    theme_void() + theme(plot.title = element_text(hjust = 0.5)) +
    ggtitle(paste0("\nBands ", paste(defbands2,collapse=", ")))
  #)
  grid.arrange(ggRad, ggRad2, ncol=2)

1.3 Polygon vector

Relevant parts of each object were interactively digitized as polygons in QGIS and used to select spectra from the image.

options(warn=-1)
nompoly <- "034.shp"
p <- vect(file.path(vecdir,nompoly))
#values(p)$Name <- c("White", "Ochre", "Ochre_800",
#                    "Red", "Red_800",
#                    "GreenTerra", "GreenTerra_800")
p

 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 3, 3  (geometries, attributes)
 extent      : 9.067334, 502.9327, 9.403755, 497.8865  (xmin, xmax, ymin, ymax)
 source      : 034.shp
 coord. ref. : Arbitrary 

#knitr::kable(p)

2. Geometric consistency

In this example, the CRS was defined in the image by editing the original hdr file, and the digitized polygons were defined in QGIS with the same CRS. We confirm the correct overlay of vector polygons and raster image (RGB display with bands 70 (598 nm), 53 (548 nm) and 19 (449 nm)).

crs(s) <- crs(p)
#ext(s) <- c(0,512,-512,0)
#plotRGB(stretch(s[[defbands2]], smin=0.0, smax=0.85), ext=c(50,400,-360,-100))
plotRGB(stretch(s[[defbands]], smin=0.0, smax=1.0), main=IQbands[defbands,])
plot(p, add=TRUE)

text(p,2,cex=0.9,font=2)

3. Object spectra

Finally, using the polygons as templates, we extract the reflectance values from the hyperspectral image and plot.

options(warn=-1)
refldat.m  <- terra::extract(s, p, fun=mean, na.rm=TRUE)
refldat.sd <- terra::extract(s, p, fun=sd, na.rm=TRUE)
refldatm <- as.data.frame(t(refldat.m[,-1]))
#colnames(refldatm) <- values(p)$Name
colnames(refldatm) <- values(p)$Name
refldatm$Band <- IQbands$Band
refldatm$Wavelength <- IQbands$Wavelength
#head(refldatm)
refldatsd <- as.data.frame(t(refldat.sd[,-1]))
colnames(refldatsd) <- values(p)$Name
#head(refldatsd)

refldat <- melt(refldatm,id.vars = c("Band", "Wavelength"))
#head(refldat)
names(refldat)[3:4] <- c("Sample", "Reflectance")
refldat$Sample <- mapvalues(refldat$Sample, from= values(p)$Name, to= values(p)$Comment)
a <- melt(refldatsd,id.vars=NULL)
#head(a)
refldat$sd <- a$value
refldat$Treatment <- mapvalues(refldat$Sample, from=values(p)$Acronym, to=values(p)$Treatment)
refldat$Pigment <- mapvalues(refldat$Sample, from=values(p)$Acronym, to=values(p)$Pigment)
#head(refldat)

#gg1 <- plothcavspect(datainRef,title="")
#Default colors
#ggplotColours <- function(n = 6, h = c(0, 360) + 15){
#  if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
# hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
#}

#micolor <- ggplotColours(n=36)
#scales:::show_col(micolor)
#https://www.nceas.ucsb.edu/sites/default/files/2020-04/colorPaletteCheatsheet.pdf
#micolor <- c("W"="cyan", "G50"="grey50", "GC"="grey20")

gg1 <- ggplot(data=refldat) +
  geom_point(aes(x=Wavelength, y=Reflectance, color=Sample),size=0.5) +
  geom_line(aes(x=Wavelength, y=Reflectance, color=Sample)) +
  geom_errorbar(aes(x=Wavelength,ymin=Reflectance-sd, ymax=Reflectance+sd, color=Sample),
                width=.001, alpha=0.5,
                position=position_dodge(0.05)) +
  ylim(c(0,1.2)) +
  xlab("Wavelength (nm)") +
  theme(legend.position = "bottom") +
  ggtitle("Grey standards")

#  gg1 + geom_dl(aes(x=Wavelength, y=Reflectance, label = Sample), method = list(dl.combine("first.points", "last.points"), cex = #0.7)) + xlim(c(300, 1100))

#+ scale_color_manual(values=micolor)

ggplotly(gg1,
         tooltip=c("Band","Wavelength",Reflectance="Reflectance","Sample"),
         height=600, width=900)

NA