A comparison of VISNIR HSI systems

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

This script extracts and compares reflectance spectra of a standard reflectance target (Labsphere WCS-MC-020) from VISNIR hyperspectral images acquired with three VISNIR hyperespectral imaging systems in HSILab premises:

• Cubert Firefleye S185 SE
• Specim IQ
  
• Specim FX10

defaultW <- getOption("warn")
options(warn=-1)
#Paths
RSpectDir       <- "/home/alobo/owncloudRSpect/RSpect"
datadirC         <- "/home/alobo/WORK/HSI_Lab/Tests_Hypercams/LlumASD-KAISER_2019-01-23/Cubert_KAISER2-3_2019-01-23"
imadirC          <- file.path(datadirC,"cue")
vecdirC          <- file.path(datadirC,"QGIS/Polygons")
datadirIQ         <- "/home/alobo/WORK/HSI_Lab/Tests_Hypercams/SpecimIQ/20220914/018"
imadirIQ          <- file.path(datadirIQ,"results")
vecdirIQ          <- file.path(datadirIQ,"QGIS/Polygons")
datadirFX10         <- "/home/alobo/WORK/HSI_Lab/Tests_Hypercams/SpecimFX10/FX10_Sant_Finx_big_2019-06-11_14-01-10/FX10/"
imadirFX10          <- file.path(datadirFX10,"capture")
vecdirFX10          <- file.path(datadirFX10,"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
load(file.path(RSpectDir,"Cubertbands.rda")) #Bands Cubert system
load(file.path(RSpectDir,"FX10bands.rda")) #Bands Cubert system
load(file.path(RSpectDir,"RefMCAA01.rda"))#Labsphere WCS-MC-020

1. Input images

1.1 Cubert Reflectance image

We use the *.cue reflectance image that was creted by Cubert’s program to export cub files that are written by Cubert software. We use the original image, with no pansharpening.

Warning: [rast] unknown extent

class       : SpatRaster 
dimensions  : 50, 50, 138  (nrow, ncol, nlyr)
resolution  : 20, 20  (x, y)
extent      : 0, 1000, -1000, 0  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
names       : 450    nm, 454     nm, 458     nm, 462     nm, 466     nm, 470     nm, ... 
min values  :   0.0043,   0.0088,   0.0027,   0.0011,   0.0236,   0.0505, ... 
max values  :   0.9410,   6.5491,   6.5532,   1.0084,   0.9868,   0.9702, ... 

1.2 Specim IQ Reflectance image

The reflectance image was calculated by the Specim IQ system itself using the average value of the White reference The original hdr file was manualy edited to include the mapinfo field as CRS information: map info = {Arbitrary, 1, 1, 0 ,512,1, 1, 0, North}

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_018.dat 
names       :    397.32,    400.20,     403.09,     405.97,     408.85,     411.74, ... 
min values  : 0.1724138, 0.1351351, 0.09803922, 0.07894737, 0.05504587, 0.03921569, ... 
max values  : 1.2142857, 1.2162162, 1.15384614, 1.11842108, 1.11009169, 1.08496737, ... 

We display all 3 RGB color composites with similar spectral bands:

1.3 Specim FX10 Reflectance image

The reflectance image was calculated by the Specim IQ system itself using the average value of the White reference The original hdr file was manualy edited to include the mapinfo field as CRS information: map info = {Arbitrary, 1, 1, 0 ,512,1, 1, 0, North}

class       : SpatRaster 
dimensions  : 1173, 1024, 448  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 1024, -149, 1024  (xmin, xmax, ymin, ymax)
coord. ref. : Arbitrary 
source      : REFL_FX10_Sant_Finx_big_2019-06-11_14-01-10.raw 
names       : REFL_~-10_1, REFL_~-10_2, REFL_~-10_3, REFL_~-10_4,     402.24, 403.55, ... 
min values  :       -0.36,  -0.3913043,       -0.28,  -0.2222222, -0.2222222,  -0.25, ... 
max values  :        2.00,   1.3333334,        1.50,   1.3750000,  1.3870968,   1.32, ... 

1.4 Display reflectance images

We display all 3 RGB color composites with similar spectral bands:

Cubert bands

	Bandnb	Band	Wavelength
138	138	Band138	998
42	42	Band042	614
26	26	Band026	550
5	5	Band005	466

IQ bands

	Band	Wavelength
203	Band203	1000.49
75	Band075	613.38
53	Band053	548.55
25	Band025	466.77

FX10 bands

	Bandnb	Band	Wavelength	fwhm
445	445	Band445	1000.29	1.41
165	165	Band165	614.55	1.35
116	116	Band116	548.91	1.33
54	54	Band054	466.64	1.32

  crs(s1) <- crs(s2) <- crs(s3) <- "" #geom_spatraster error otherwise
  ggRad <- ggplot() + 
  geom_spatraster_rgb(data = stretch(subset(s1, subset=defbandsC), minq=0.0, maxq=1)) +
    coord_fixed(ratio = 1) +
    theme_void() + theme(plot.title = element_text(hjust = 0.5)) +
    ggtitle(paste0("Cubert S185 SE \nBands ", paste(Cubertbands$Wavelength[defbandsC],collapse=", ")))
  ggRad2 <- ggplot() + 
    geom_spatraster_rgb(data = stretch(subset(s2, subset=defbandsIQ), minq=0.0, maxq=1)) +
    coord_fixed(ratio = 1) +
    theme_void() + theme(plot.title = element_text(hjust = 0.5)) +
     ggtitle(paste0("Specim IQ \nBands ", paste(IQbands$Wavelength[defbandsIQ],collapse=", ")))
   ggRad3 <- ggplot() + 
    geom_spatraster_rgb(data = stretch(subset(s3, subset=defbandsFX10), minq=0.0, maxq=0.95)) +
    coord_fixed(ratio = 1) +
    theme_void() + theme(plot.title = element_text(hjust = 0.5)) +
     ggtitle(paste0("Specim FX10 \nBands ", paste(FX10bands$Wavelength[defbandsFX10],collapse=", ")))

  grid.arrange(ggRad, ggRad2, ggRad3, ncol=3)

2. Polygon vectors

Relevant parts of each reference target were interactively digitized as polygons in QGIS and used to select spectra from the images.

options(warn=-1)
nompolyC <- "KAISER2-3.shp"
p1 <- vect(file.path(vecdirC,nompolyC))
#values(p)$Name <- c("White", "Ochre", "Ochre_800",
#                    "Red", "Red_800",
#                    "GreenTerra", "GreenTerra_800")
values(p1)$Name <- c("WCS-MC-020", "White")
nompolyIQ <- "018.shp"
p2 <- vect(file.path(vecdirIQ,nompolyIQ))
#values(p2)
nompolyFX10 <- "FX10_Sant_Finx_big_2019-06-11_14-01-10.shp"
p3 <- vect(file.path(vecdirFX10,nompolyFX10))

crs(s1) <- crs(p1)
plotRGB(stretch(s1[[defbandsC]], smin=0.0, smax=1),
        main=Cubertbands[defbandsC,])
#lines(p1, xlim=c(50,200))
lines(p1)

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

crs(s2) <- crs(p2)
plotRGB(stretch(s2[[defbandsIQ]], smin=0.0, smax=1),
        main=IQbands[defbandsIQ,])
#lines(p2, xlim=c(50,200))
lines(p2)

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

crs(s3) <- crs(p3)
plotRGB(stretch(s3[[defbandsFX10]], smin=0.0, smax=1),
        main=IQbands[defbandsFX10,])
#lines(p2, xlim=c(50,200))
lines(p3)

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

3. Object spectra

Using the polygons as templates, we extract the reflectance values from the hyperspectral images.

3.1 Cubert

options(warn=-1)
refldat1.m  <- terra::extract(s1, p1, fun=mean, na.rm=TRUE)
refldat1.sd <- terra::extract(s1, p1, fun=sd, na.rm=TRUE)

refldatm1 <- as.data.frame(t(refldat1.m[,-1]))
colnames(refldatm1) <- values(p1)$Name
refldatm1$Band <- Cubertbands$Band
refldatm1$Wavelength <- Cubertbands$Wavelength
#head(refldatm)
refldatsd1 <- as.data.frame(t(refldat1.sd[,-1]))
colnames(refldatsd1) <- values(p1)$Name
#head(refldatsd)

refldat1 <- melt(refldatm1,id.vars = c("Band", "Wavelength"))
#head(refldat)
names(refldat1)[3:4] <- c("Sample", "Reflectance")
a <- melt(refldatsd1,id.vars=NULL)
#head(a)
refldat1$sd <- a$value
refldat1$Source = "Cubert S185SE"
#head(refldat1)

3.2 Specim IQ

options(warn=-1)
refldat2.m  <- terra::extract(s2, p2, fun=mean, na.rm=TRUE)
refldat2.sd <- terra::extract(s2, p2, fun=sd, na.rm=TRUE)

refldatm2 <- as.data.frame(t(refldat2.m[,-1]))
colnames(refldatm2) <- values(p2)$Name
refldatm2$Band <- IQbands$Band
refldatm2$Wavelength <- IQbands$Wavelength
#head(refldatm)
refldatsd2 <- as.data.frame(t(refldat2.sd[,-1]))
colnames(refldatsd2) <- values(p2)$Name
#head(refldatsd)

refldat2 <- melt(refldatm2,id.vars = c("Band", "Wavelength"))
#head(refldat2)
names(refldat2)[3:4] <- c("Sample", "Reflectance")
a <- melt(refldatsd2,id.vars=NULL)
#head(a)
refldat2$sd <- a$value
refldat2$Source = "Specim IQ"
#head(refldat2)

3.3 Specim FX10

options(warn=-1)
refldat3.m  <- terra::extract(s3, p3, fun=mean, na.rm=TRUE)
refldat3.sd <- terra::extract(s3, p3, fun=sd, na.rm=TRUE)

refldatm3 <- as.data.frame(t(refldat3.m[,-1]))
colnames(refldatm3) <- values(p3)$Name
refldatm3$Band <- FX10bands$Band
refldatm3$Wavelength <- FX10bands$Wavelength
#head(refldatm)
refldatsd3 <- as.data.frame(t(refldat3.sd[,-1]))
colnames(refldatsd3) <- values(p3)$Name
#head(refldatsd)

refldat3 <- melt(refldatm3,id.vars = c("Band", "Wavelength"))
#head(refldat2)
names(refldat3)[3:4] <- c("Sample", "Reflectance")
a <- melt(refldatsd3,id.vars=NULL)
#head(a)
refldat3$sd <- a$value
refldat3$Source = "Specim FX10"
#head(refldat2)

3.4 Combine

We combine all 3 datasets and organize the certified spectra of the Labsphere WCS MC-020 target in an equivalent dataset

refldat <- rbind(refldat1, refldat2, refldat3)
refldat$Reflectance[refldat$Reflectance>5] <- NA
a <- data.frame( Band=NA, Wavelength=RefMCAA01$Wavelength, Sample="WCS-MC-020", 
                 Reflectance=RefMCAA01$Reflectance, sd=0.0, Source="Certified")

4. Plot reflectance spectra

We plot the spectra, including the certified spectra of the Labsphere WCS MC-020 target:

options(warn=-1)
#Custom colors from Default colors
#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)
#scales:::show_col(micolor)
#https://www.nceas.ucsb.edu/sites/default/files/2020-04/colorPaletteCheatsheet.pdf
micolor <- c("Certified"="black", "Cubert S185SE"= micolor[4], 
             "Specim FX10"=micolor[2], "Specim IQ" = micolor[3])

gg1 <- ggplot() +
  #geom_line(data=a, aes(x=Wavelength, y=Reflectance), color="black", size=0.5) +
  geom_point(aes(x=Wavelength, y=Reflectance, color=Source),size=0.5) +
  geom_line(aes(x=Wavelength, y=Reflectance, color=Source)) +
  geom_errorbar(aes(x=Wavelength,ymin=Reflectance-sd, ymax=Reflectance+sd, color=Source),
                width=.001, alpha=0.5,
                position=position_dodge(0.05)) +
  xlab("Wavelength (nm)") +
  xlim(c(400, 1000)) +
  scale_color_manual(values=micolor) +
  theme(legend.position = c(0.85, 0.22)) 

gg2 <- gg1 %+% rbind(refldat[refldat$Sample=="WCS-MC-020",],a) +
  ggtitle("Comparison of HSI systems on Labsphere WCS-MC-020 target")

gg3 <- gg1 %+% refldat[refldat$Sample=="White",] +
  ggtitle("Comparison of HSI systems on White target")

#gg1 + facet_wrap(~Sample,nrow=2) 

ggplotly(gg2 + facet_wrap(~Sample,nrow=2) ,
         tooltip=c("Band","Wavelength",Reflectance="Reflectance","Source"),
         height=600, width=1200)

ggplotly(gg3 + facet_wrap(~Sample,nrow=2) ,
         tooltip=c("Band","Wavelength",Reflectance="Reflectance","Source"),
         height=600, width=1200)