Stability of SpectroRasp spectra

• Agustin.Lobo@geo3bcn.csic.es
• 20221222, 20230103, 20230109, 20230206

Click on “Code” symbols to unfold code (folded by default).

require(ggplot2)
require(grid)
require(gridExtra)
require(gganimate)
require(plotly)
require(lubridate)
require(colorRamps)
require(lmodel2)
require(broom)
require(pander)
source("selSessionRefs.R")
#Certified
load("/home/alobo/owncloudRSpect/RSpect/RefMCAA01.rda")
#load("C:/Users/Specim/ownCloud/RSpect/RefMCAA01.rda")

1. Stability of White Reference spectra

1.1 Stability of White Reference Radiance spectra (Data on cont. measurements 20221220)

We assess the stability of HSILab SpectroPi readings using the Spectralon White reference Method: continuous (every minute) radiance readings using time.measure.py with settings:

• IT STS: 126 ms
• IT NIRQUEST: 30 ms
• Scans: 20

Radiance readings are not stable in the long term, but reading the white reference prior to each sample (as we do in our protocol) successfully results on stable and reproducible reflectance measurements.

fdir <- "../USBRaspberryBackup_20221219/time_measures2022"
lista <- data.frame(fname=list.files(fdir, patt=glob2rx("*.txt")), fdir=fdir)
Radcont <- data.frame(V1=0, V2=0, Radianceant=0, RadianceDif0=0, RadianceDif=0, sample=NA)
for (i in 1:nrow(lista)){
  if(i >1) radant <- a$V2
  a <- read.csv(file.path(lista[i,2], lista[i,1]),sep=" ", header=FALSE)
  if(i==1) radini <- radant <- a$V2
  a$Radianceant <- radant
  a$RadianceDif0 <- a$V2-radini
  a$RadianceDif <- a$V2-radant
  a$sample <- i
  Radcont <- rbind(Radcont,a)
}
colnames(Radcont) <- c("Wavelength", "Radiance", "Radianceant", "RadianceDif0","RadianceDif","Sample")
Radcont <- Radcont[-1,]
# dim(Radcont)
# head(Radcont)

#Display Selected readings:
sel <- Radcont[Radcont$Sample %in% c(2,25,50,75,100,135),]
sel$Reflectance0 <- sel$Radiance/Radcont$Radiance[Radcont$Sample==1]
sel$Reflectanceant <- sel$Radiance/sel$Radianceant
#table(sel$Sample)
samplenames <- paste0("t=",unique(sel$Sample))
names(samplenames) <- unique(sel$Sample)

#Raw Radiance data
ggRadcontsel <- ggplot(data=sel) +
  geom_point(aes(x=Wavelength, y=Radiance, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  geom_line (aes(x=Wavelength, y=Radiance, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  xlim(c(400, 2600))  + 
  #facet_wrap(~Sample, labeller=as_labeller(samplenames)) + 
  theme(legend.position="none") +
  #ggtitle("Stability of Radiance (White Reference)", subtitle = "Raw data")
  ggtitle("Radiance", subtitle = "Raw data")
#ggRadcontsel
ggsave(ggRadcontsel, filename="ggRadcontsel.jpeg",units="px", width=2100, height=2100)

#Radiance data (1st reading subtracted)
ggRadDif0contsel <- ggplot(data=sel) +
  geom_point(aes(x=Wavelength, y=RadianceDif0, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  geom_line (aes(x=Wavelength, y=RadianceDif0, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  xlim(c(400, 2600)) + ylim(c(-200,600)) + 
  #facet_wrap(~Sample, labeller=as_labeller(samplenames)) +
  ylab("Radiance - Radiance at t=0") +  
  theme(legend.position="none") +
  #ggtitle("Stability of Radiance (White Reference)", subtitle = "initial Radiance subtracted")
  ggtitle("Radiance", subtitle = "initial Radiance subtracted")
#ggRadDif0contsel

#grid.arrange(ggRadcontsel, ggRadDif0contsel, ncol=2,
#                           top=textGrob("Stability of Readings (White Reference)"))
ggsave(ggRadDif0contsel, filename="ggRadDif0contsel.jpeg",units="px", width=2100, height=2100)

#Display all radiance readings (raw data)
ggRadcont <- ggplot(data=Radcont,aes(x=Wavelength, y=Radiance, group=Sample,color=Sample)) +
  geom_point(alpha=1, size=0.5) +
  geom_line (alpha=1, linewidth=0.5) +
  #geom_line(data=Radcont[Radcont$Sample==2,]) +
  xlim(c(400, 2600))  +
  ggtitle("Stability of Radiance (White Reference)", subtitle = "Raw data")
#Visualizing all readings (samples) at once is not useful
# ggRadcont
# #ggplotly(ggRadcont) takes too long
# Rather use animation:

#not run in notebook:
#ggRadcont.anim <- ggRadcont + transition_time(Sample) +
#  labs(title = "Time: {frame_time}")
#animate(ggRadcont.anim)
#anim_save("ggWRadcont.anim.gif")

#Display all radiance readings (1st reading subtracted)
ggRadDifcont <- ggplot(data=Radcont, aes(x=Wavelength, y=RadianceDif0, group=Sample,color=Sample)) +
  #geom_point(alpha=1, size=0.5) +
  #geom_line (alpha=1, linewidth=0.5) +
  geom_line(data=Radcont[Radcont$Sample==2,]) +
  xlim(c(400, 2600))  +
  ggtitle("Stability of Radiance (White Reference)", subtitle = "initial Radiance subtracted")
#ggRadDifcont

#not run in notebook:
#ggRadDifcont_anim1 <- ggRadDifcont + transition_time(Sample) +
#  labs(title = "Time: {frame_time}") 
#animate(ggRadDifcont_anim1)
#anim_save("ggWRadDifcont_anim1.gif")

Raw Data	Initial reading subtracted
animate	animate

1.2 Stability of White reference Reflectance spectra (Data on cont. measurements 20221220)

The variablity of Reflectance readings is compensated by the ratio to a radiance reading that is close in time (eg. immediately prior) to that of the sample, as it is customary in lab and field protocols.

#Display Selected readings:
fdir <- "../USBRaspberryBackup_20221219/time_measures2022"
sel <- Radcont[Radcont$Sample %in% c(2,25,50,75,100,135),]
sel$Reflectance0 <- sel$Radiance/Radcont$Radiance[Radcont$Sample==1]
sel$Reflectanceant <- sel$Radiance/sel$Radianceant
#table(sel$Sample)
samplenames <- paste0("t=",unique(sel$Sample))
names(samplenames) <- unique(sel$Sample)

#Raw Radiance data
ggRadcontsel <- ggplot(data=sel) +
  geom_point(aes(x=Wavelength, y=Radiance, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  geom_line (aes(x=Wavelength, y=Radiance, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  xlim(c(400, 2600))  + facet_wrap(~Sample, labeller=as_labeller(samplenames)) + theme(legend.position="none") +
  #ggtitle("Stability of Radiance (White Reference)", subtitle = "Raw data")
  ggtitle("Radiance", subtitle = "Raw data")
#ggRadcontsel

#Radiance data (1st reading subtracted)
ggRadDif0contsel <- ggplot(data=sel) +
  geom_point(aes(x=Wavelength, y=RadianceDif0, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  geom_line (aes(x=Wavelength, y=RadianceDif0, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  xlim(c(400, 2600)) + ylim(c(-200,600)) + facet_wrap(~Sample, labeller=as_labeller(samplenames)) +
  ylab("Radiance - Radiance at t=0") +  theme(legend.position="none") +
  #ggtitle("Stability of Radiance (White Reference)", subtitle = "initial Radiance subtracted")
  ggtitle("Radiance", subtitle = "initial Radiance subtracted")
#ggRadDif0contsel

#Reflectance data (basis: 1st radiance reading)
ggRefl0contsel <- ggplot(data=sel) +
  geom_point(aes(x=Wavelength, y=Reflectance0, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  geom_line (aes(x=Wavelength, y=Reflectance0, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  xlim(c(400, 2600)) + ylim(c(0.9, 1.1)) + facet_wrap(~Sample, labeller=as_labeller(samplenames)) + 
  ylab("Radiance/Radiance at t=0") +  theme(legend.position="none") +
  #ggtitle("Stability of Reflectance (White Reference)", subtitle = "basis: initial Radiance")
  ggtitle("Reflectance", subtitle = "basis: initial Radiance")
#ggRefl0contsel

#Reflectance data (basis: previous radiance reading)
ggReflantcontsel <- ggplot(data=sel) +
  geom_point(aes(x=Wavelength, y=Reflectanceant, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  geom_line (aes(x=Wavelength, y=Reflectanceant, group=Sample,color=Sample),alpha=0.5, size=0.5) +
  xlim(c(400, 2600)) +  ylim(c(0.9, 1.1)) + facet_wrap(~Sample, labeller=as_labeller(samplenames)) +
  ylab("Radiance/Radiance at t=t-1") +  theme(legend.position="none") +
  #ggtitle("Stability of Reflectance (White Reference)", subtitle = "basis: previous Radiance")
  ggtitle("Reflectance", subtitle = "basis: previous Radiance")
#ggReflantcontsel

Wstability <- grid.arrange(ggRadcontsel, ggRadDif0contsel, ggRefl0contsel,ggReflantcontsel,
                           top=textGrob("Stability of Readings (White Reference)"))

#print(Wstability)
#not run in nb
#ggsave(Wstability, file="Wstability.pdf", width=12, height=8)

Selecting a given wavelength, we emphasize the drift of radiance readings and the stabilization by using the ratio to the previous (within 1’) reading: reflectance variability becomes bounded to +/- 0.0001

#Variability of Reflectance @1418 nm 
Radcont1418 <- Radcont[round(Radcont$Wavelength)==1419 ,]
Radcont1418$Reflectance0 <- Radcont1418$Radiance/Radcont1418$Radiance[Radcont1418$Sample==1]
Radcont1418$Reflectanceant <- Radcont1418$Radiance/Radcont1418$Radianceant

#Radiance @1418 nm
ggRad1418 <- ggplot(data=Radcont1418) +
  geom_point(aes(x=as.numeric(Sample), y=Radiance)) +
  xlab("Time (min)") + ylab("Radiance@1418.5820 nm")+ 
  ggtitle("Radiance")
#ggRad1418

#Reflectance @1418 nm  (basis: 1st radiance reading)
ggRefl01418 <- ggplot(data=Radcont1418) +
  geom_point(aes(x=as.numeric(Sample), y=Reflectance0)) +
  xlab("Time (min)") + ylab("Reflectance@1418.5820 nm")+ 
  ylim(c(0.999, 1.005)) +
  ggtitle("", subtitle = "basis: initial Radiance")
#ggRefl01418

#Reflectance @1418 nm  (basis: previous radiance reading)
ggReflant1418 <- ggplot(data=Radcont1418) +
  geom_point(aes(x=as.numeric(Sample), y=Reflectanceant)) +
  xlab("Time (min)") + ylab("Reflectance@1418.5820 nm")+ 
  ylim(c(0.999, 1.005)) +
  ggtitle("", subtitle = "basis: previous Radiance")
#ggReflant1418

#Wstability1418 <- grid.arrange(ggRad1418,ggRefl01418,ggReflant1418)
#ggsave(Wstability1418, file="Wstability1418.pdf", width=10, height=14, top=textGrob("Stability of Readings (White Reference)"))
Wstability1418 <- grid.arrange(ggRefl01418,ggReflant1418)

2. Stability of Standard Reference spectra (20230102 data)

Initial and final (138’ later) Reflectance measurements of the Standard reference are slightly shifted, despite having measured the white reference immediately before the target. Integration times (IT) were set to those calculated as optima by SpR: 46 ms (STS) and 15 ms (NIRQUEST).

fdir <- "../time_measures20230102"
lista <- data.frame(fname=list.files(fdir, patt=glob2rx("*.txt")), fdir=fdir)[-(1:121),] #exclude auto-capture readings
#lista <- lista[-1,]

#Inital and final (138' later) Reflectance measurements with IT 46 ms (STS) and 15 ms (NIRQUEST)
listasel <- lista[1:2,]
selRefs <- selSessionRefs(selflist=listasel)
selRefs$sample[selRefs$sample=="Ref_1"] <- "Initial" 
selRefs$sample[selRefs$sample=="Ref_2"] <- "Final" 
ggRefRefl <- ggplot(data=selRefs,aes(x=Wavelength, y=Reflectance)) +
  geom_point(alpha=1, size=0.5,aes(group=sample,color=sample)) +
  geom_line (alpha=1, linewidth=0.5, aes(group=sample, color=sample)) +
  #geom_line(data=Radcont[Radcont$Sample==2,]) +
  #geom_line(data=RefMCAA01, aes(x=Wavelength, y=Reflectance),col="grey") +
  geom_vline(aes(xintercept=1418.5820),linetype=4) +
  xlim(c(400, 2600)) +
  theme(legend.position="bottom") +
  ggtitle("Reflectance (Standard Target)", subtitle = "Initial and final readings (IT: 46 ms, 15 ms)")
#ggplotly(ggRefRefl)
ggplotly(ggRefRefl)  %>%
  layout(legend = list(orientation = "h", x = 0.7, y = 0.1))

Initial and Final Radiance and Reflectance spectra are linearly related with slopes close to 1. Relatively large intercept values indicate again a drift in Radiance, which is reduced (but not eliminated) in Reflectance measurements: 0.006 for the STS-IT (400-900 nm) and -0.07 for the NIRQUEST (950-2600 nm).

d <- selRefs[,c(1:6)]
ds <- plyr::ddply(d, c("sample","Wavelength"), summarise, 
                  Radiance = (Radiance-Dark),
                  WRadiance = WhiteRad,
                  Reflectance = Reflectance)

dsw <- data.frame(Wavelength=ds$Wavelength[ds$sample=="Initial"],
                  InitialRefl=ds$Reflectance[ds$sample=="Initial"], FinalRefl=ds$Reflectance[ds$sample=="Final"],
                  Initial=ds$Radiance[ds$sample=="Initial"], Final=ds$Radiance[ds$sample=="Final"],
                  WInitial=ds$WRadiance[ds$sample=="Initial"], WFinal=ds$WRadiance[ds$sample=="Final"])
dsw$Instrument <- "STS-IT"
dsw$Instrument[dsw$Wavelength>850] <- "NIRQUEST"

lm2Rad1 <- lmodel2(Final~Initial, data=dsw[dsw$Instrument=="STS-IT",])
lm2Rad2 <- lmodel2(Final~Initial, data=dsw[dsw$Instrument=="NIRQUEST",])
#lm2Rad1$regression.results[2,]
#lm2Rad2$regression.results[2,]

lm2W1 <- lmodel2(WFinal~WInitial, data=dsw[dsw$Instrument=="STS-IT",])
lm2W2 <- lmodel2(WFinal~WInitial, data=dsw[dsw$Instrument=="NIRQUEST",])
#lm2W1$regression.results[2,]
#lm2W2$regression.results[2,]

lm2Refl1 <- lmodel2(FinalRefl~InitialRefl, data=dsw[dsw$Instrument=="STS-IT",])
lm2Refl2 <- lmodel2(FinalRefl~InitialRefl, data=dsw[dsw$Instrument=="NIRQUEST",])
#lm2Refl1$regression.results[2,]
#lm2Refl2$regression.results[2,]

a <- rbind(
  lm2W1$regression.results[2,2:3],
  lm2W2$regression.results[2,2:3],
  lm2Rad1$regression.results[2,2:3],
  lm2Rad2$regression.results[2,2:3],
  lm2Refl1$regression.results[2,2:3],
  lm2Refl2$regression.results[2,2:3])
a$Instrument <- rep(c("STS-IT", "NIRQUEST"),3)
a$Target <- c(rep("W",2), rep("Std",4))
a$Magnitude <- c(rep("Radiance",4), rep("Reflectance",2))
a <- a[,c(4,5,3,1,2)]
rownames(a) <- NULL


ggWRad <- ggplot(data=dsw, aes(x=WInitial, y=WFinal,color=Instrument)) +
  geom_point() +
  geom_abline(aes(intercept=lm2W1$regression.results[2,2], slope=lm2W1$regression.results[2,3])) +
  geom_abline(aes(intercept=lm2W2$regression.results[2,2], slope=lm2W2$regression.results[2,3])) +
  geom_abline(aes(intercept=0, slope=1), linetype=4) +
  xlab("Radiance (Initial measurement)") +  ylab("Radiance (Final measurement)") + 
  theme(legend.position = c(.85,.15), aspect.ratio=1) + 
  ggtitle("White Reference") 
ggRad <- ggplot(data=dsw, aes(x=Initial, y=Final,color=Instrument)) +
  geom_point() +
  geom_abline(aes(intercept=lm2Rad1$regression.results[2,2], slope=lm2Rad1$regression.results[2,3])) +
  geom_abline(aes(intercept=lm2Rad2$regression.results[2,2], slope=lm2Rad2$regression.results[2,3])) +
  geom_abline(aes(intercept=0, slope=1), linetype=4) +
  xlab("Radiance (Initial measurement)") +  ylab("Radiance (Final measurement)") + 
  theme(legend.position = c(.85,.15), aspect.ratio=1) +
  ggtitle("Standard Reference")
ggRefl <- ggplot(data=dsw, aes(x=InitialRefl, y=FinalRefl,color=Instrument)) +
  geom_point() +
  geom_abline(aes(intercept=lm2Refl1$regression.results[2,2], slope=lm2Refl1$regression.results[2,3])) +
  geom_abline(aes(intercept=lm2Refl2$regression.results[2,2], slope=lm2Refl2$regression.results[2,3])) +
  geom_abline(aes(intercept=0, slope=1), linetype=4)+
  xlab("Reflectance (Initial measurement)") +  ylab("Reflectance (Final measurement)") + 
  theme(legend.position = c(.85,.15), aspect.ratio=1) +
  ggtitle("Standard Reference")
grid.arrange(ggWRad, ggRad,  nrow=1)

pander::pander(arrange(a[1:4,]))#dplyr::arrange returns the df without rownames

------------------------------------------------------
 Target   Magnitude   Instrument   Intercept   Slope  
-------- ----------- ------------ ----------- --------
   W      Radiance      STS-IT       37.3      1.033  

   W      Radiance     NIRQUEST     -64.54     0.9872 

  Std     Radiance      STS-IT       37.26     1.018  

  Std     Radiance     NIRQUEST     -63.48     0.9601 
------------------------------------------------------

ggRefl

pander::pander(arrange(a[5:6,]))

--------------------------------------------------------
 Target    Magnitude    Instrument   Intercept   Slope  
-------- ------------- ------------ ----------- --------
  Std     Reflectance     STS-IT     0.006145    0.9826 

  Std     Reflectance    NIRQUEST    -0.07384    1.067  
--------------------------------------------------------