Specim IQ test with art pigments (20180918)

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

This script extracts reflectance spectra of a set of art pigments from a VISNIR hyperspectral image acquired with a demo unit of the Specim IQ system in HSILab premises on 20180918:

• Specim IQ white reference card.
• Ocre (Ochre)
• Ocre 800 (Ochre heated to 800ºC)
• Vermell (Iron Oxide Red)
• Red 800 (Iron Oxide Red heated to 800ºC)
• Terra (green earth)
• Terra Verda 800 (green earth heated to 800ºC)

Pigment samples were provided by jibanez@geo3bcn.csic.es

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/2018-10-18_008"
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_2018-10-18_008.png"))
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_2018-10-18_008.png 
names       : RGBBACK~8_008_1, RGBBACK~8_008_2, RGBBACK~8_008_3, RGBBACK~8_008_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. :  
source      : REFLECTANCE_2018-10-18_008.dat 
names       : 397.32, 400.20, 403.09, 405.97, 408.85, 411.74, ... 

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 <- "008.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  : 7, 7  (geometries, attributes)
 extent      : 72.29907, 382.5495, -352.3252, -180.9495  (xmin, xmax, ymin, ymax)
 source      : 008.shp
 coord. ref. : ETRS89 / UTM zone 31N (EPSG:25831) 

#knitr::kable(p)

2. Make consistent geometry

In case the CRS is not defined in the image (and this is the unfortunate case of Specim hdr files), R and QGIS assume, by default, different geometry for the same tiff. As a consequence, polygons digitized in QGIS do not overlay the raster image in R. The code must account for this fact, making a consistent geometry for raster and vector layer, which we confirm in the 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=0.4), ext=c(50,200,-400,-150),
        main=IQbands[defbands,])
lines(p, xlim=c(50,200))

text(p,5,cex=0.9,font=2)
text(p,5,cex=0.9,font=1,col=0)

#lines(p,ext=c(50,400,-360,-100)) # "ext" is not a graphical parameter
#lines(ext(c(50,400,-360,-100)))

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

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"="grey70", 
             "Och"="coral4",
             "Och_800"="chocolate4",
             "R"="red",
             "R_800"="red4",
             "G"="olivedrab3",
             "G_800"="olivedrab")

gg1 <- ggplot(data=refldat) +
  geom_point(aes(x=Wavelength, y=Reflectance, color=Sample, shape=Treatment),size=0.5) +
  geom_line(aes(x=Wavelength, y=Reflectance, color=Sample, linetype=Treatment)) +
  geom_errorbar(aes(x=Wavelength,ymin=Reflectance-sd, ymax=Reflectance+sd, color=Sample),
                width=.001, alpha=0.5,
                position=position_dodge(0.05)) +
  xlab("Wavelength (nm)") +
  theme(legend.position = "bottom") +
  ggtitle("Art pigments")

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

#gg1 + xlim(c(450, 1000))  + 
#  facet_wrap(~State) + geom_dl(aes(x=Wavelength, y=Reflectance, label = Sample), method = list(dl.combine("first.points", "last.points"), cex = 0.7)) 
#gg2 <- gg1 + 

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

NA