Specim IQ test with lava sample and standard

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

This script extracts reflectance spectra of 3 objects from a VISNIR hyperspectral image acquired by the Specim IQ system in HSILab premises:

• Grey Reference Zenith Lite SG 3153 50%
• Labsphere WCS-MC-020 standard.
• DIPV-23 (IED180023) sample (“famous” lava flow body entering into the lake).
  • Crater Lake (Mt. Kikwood) Deception Is. (Antarctica) 201802dd).

As the DIPV-23 sample is very dark, we used a 50% grey reference instead of a white to adjust integration time to more appropriate values while avoiding saturation (The WCS standard is quite bright, thus no big advantage in this example).

Code

defaultW <- getOption("warn")
options(warn=-1)
#Paths
RSpectDir       <- "/home/alobo/owncloudRSpect/RSpect"
datadir         <- "/home/alobo/WORK/HSI_Lab/Tests_Hypercams/SpecimIQ/SpecimIQTest_2018-09-18/IQ_DIPV23_20180918"
imadir          <- file.path(datadir,"results")
vecdir          <- "/home/alobo/WORK/HSI_Lab/Tests_Hypercams/SpecimIQ/SpecimIQTest_2018-09-18/IQ_DIPV23_20180918/QGIS/Polygons"
#libraries
library(reshape2, quietly = TRUE)
library(plyr, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(plotly, quietly = TRUE)
library(rgdal, quietly = TRUE)
library(raster, quietly = TRUE)
#Required information
load(file.path(RSpectDir,"RefMCAA01.rda"))#Labsphere WCS-MC-020
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.

Code

nompan <- "RGBBACKGROUND_2018-10-18_004.png"
r <- brick(file.path(imadir,nompan))
r

class      : RasterBrick 
dimensions : 645, 645, 416025, 4  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 645, 0, 645  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : RGBBACKGROUND_2018-10-18_004.png 
names      : RGBBACKGROUND_2018.10.18_004.1, RGBBACKGROUND_2018.10.18_004.2, RGBBACKGROUND_2018.10.18_004.3, RGBBACKGROUND_2018.10.18_004.4 
min values :                              5,                              8,                              0,                              0 
max values :                            255,                            255,                            243,                            255 

Code

plotRGB(r)

Code

#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 Grey reference (which is interactively selected by the user during the acquisition process). The code takes the 50% reflectance of the Grey into account.

Code

nomref <- list.files(imadir,patt='dat')[1]
#nompan <- unlist(strsplit(nompan,"[.]"))[1]
s <- brick(file.path(imadir,nomref))
#s <- s * wfactor/10000
wfactor   <- 0.50 #  grey was used as white reference
s <- s * wfactor  #IQ reflectance in floating in range 0,1 
defbands <- c(70,53,19)
s

class      : RasterBrick 
dimensions : 512, 512, 262144, 204  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 512, 0, 512  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      :     X397.32,     X400.20,     X403.09,     X405.97,     X408.85,     X411.74,     X414.63,     X417.52,     X420.40,     X423.29,     X426.19,     X429.08,     X431.97,     X434.87,     X437.76, ... 
min values : 0.071428575, 0.062500000, 0.050000001, 0.059999999, 0.058823530, 0.056818184, 0.053571429, 0.050000001, 0.042168673, 0.036842104, 0.037735850, 0.030973451, 0.032258064, 0.026315790, 0.027397260, ... 
max values :   0.8461539,   0.8437500,   0.8571429,   0.8653846,   0.8333333,   0.8295454,   0.8214286,   0.7857143,   0.7771084,   0.7842105,   0.8113208,   0.8640351,   0.9032258,   0.9172933,   0.9000000, ... 

1.3 Polygon vector

Significant parts of each object were interactively digitized as polygons in QGIS and used as input here.

Code

options(warn=-1)
nompoly <- list.files(vecdir,patt='shp')
nompoly <- unlist(strsplit(nompoly,"[.]"))[1]
p <- readOGR(dsn=vecdir, layer=nompoly, verbose=FALSE)
#Use better names for samples:
p@data$Name <- c("Grey", "WCS-MC-020", "DIPV23")
spplot(p,zcol=1,colorkey=FALSE)

Code

p

class       : SpatialPolygonsDataFrame 
features    : 3 
extent      : 92.39626, 494.6542, -428.7103, -66.33271  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 3
names       : id,       Name,  Descript 
min values  :  1,     DIPV23,    DIPV23 
max values  :  3, WCS-MC-020, Reference 

2. Make consistent Geometry

In case CRS is not defined (as it is the unfortunate case of Specim hdr files), R and QGIS assume different geometry for the same tiff, thus digitized polygons do not overlay the raster image. The code must account for this fact, making a consistent geometry for raster and vector layer.

Code

#projection(r) <- projection(s) <- projection(p) #no pan image r
projection(s) <- projection(p)
extent(s) <- c(0,512,-512,0)
#plotRGB(s[[defbands]], stretch="lin") #min max different for each band
#plotRGB(stretch(s[[defbands]], minq=0.5, maxq=0.95))  #min max different for each band
plotRGB(stretch(s[[defbands]], smin=0.0, smax=1))
lines(p, col="red")

3. Extract spectra corresponding to polygons

From the next graphics, we can conclude:

1. The standard deviation of the Grey reference is quite high, indicating the occurrence of heteregeneous illumination that is not accounted for by the Specim IQ software, which uses the global average of the reference. Care must be taken to accomplish a more uniform illumination during acquisition.
2. Note that the band width (fwhm) of the IQ instrument implies that narrow absorption features in the certified spectrum appear smoother in the IQ spectra.

Code

options(warn=-1)
refldat.m  <- extract(s, p, fun=mean, na.rm=TRUE, df=TRUE)
refldat.sd <- extract(s, p, fun=sd, na.rm=TRUE, df=TRUE)

refldatm <- as.data.frame(t(refldat.m[,-1]))
colnames(refldatm) <- p@data$Name
refldatm$Band <- IQbands$Band
refldatm$Wavelength <- IQbands$Wavelength
#head(refldatm)
refldatsd <- as.data.frame(t(refldat.sd[,-1]))
colnames(refldatsd) <- p@data$Name
#head(refldatsd)

refldat <- melt(refldatm,id.vars = c("Band", "Wavelength"))
#head(refldat)
names(refldat)[3:4] <- c("Sample", "Reflectance")

a <- melt(refldatsd,id.vars=NULL)
#head(a)
refldat$sd <- a$value
#head(refldat)

a <- data.frame(Band=NA, Wavelength=RefMCAA01$Wavelength,
                Sample="Certification", 
                Reflectance=RefMCAA01$Reflectance, sd=0.0)
refldat <- rbind(refldat,a)

#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=4)
micolor <- c("Grey"="grey50", "WCS-MC-020"=micolor[4],
             "DIPV23"=micolor[3], "Certification"="black")
gg1 <- ggplot(data=refldat) +
  geom_point(aes(x=Wavelength, y=Reflectance, color=Sample),size=1) +
  geom_line(aes(x=Wavelength, y=Reflectance, color=Sample)) +
  geom_errorbar(aes(x=Wavelength,ymin=Reflectance-sd, ymax=Reflectance+sd, color=Sample),
                width=.01,
                position=position_dodge(0.05)) +
  xlab("Wavelength (nm)") +
  theme(legend.position = c(0.9, 0.30)) +
  scale_color_manual(values=micolor) +
  ggtitle("Specim IQ Test 20180918_DIPV23")
gg2 <- gg1 + geom_rect(aes(xmin = 600, xmax = 700, ymin = 0.55, ymax =1.0),col="red",alpha=0) + xlim(c(400, 1000))
print(gg2)

Code

gg1 + xlim(c(600, 700)) + ylim(c(0.55,1.0))

Code

ggplotly(gg2,
         tooltip=c("Band","Wavelength",Reflectance="Reflectance","Sample"))