Extracting reflectance data from SpectraSuite JCAMP files in R

[This article was first published on PlanetFlux, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I’ve been using an OceanOptics USB4000 spectrometer in research on biomass accumulation and climate in South Africa.

This post presents a R function that will read in JCAMP files from the SpectraSuite software and output a vector of reflectance values.  SpectraSuite can convert binary output files to JCAMP files if they weren’t saved as JCAMP.



spectracompile <- function(x,justwavelength=F,whitestandard="",darkstandard="",clip=T,range=c(500,950)) {
## Read in a .jdx file from SpectraSuite and calculate % Reflectance
## x is the filename
## justwavelengths is a logical indicating if you want to just return the vector of wavelengths
## whitestandard is an optional vector of white reference values to use if you want to override the ones in the file
## darkstandard is an optional vector of dark reference values to use if you want to override the ones in the file
## clip is a logical indicating if you want to clip the output to 0-100% reflectance. Occasionally values in noisy areas are very large
## range is the range of wavelengths you want to keep. This is useful to discard the noisy edges of the spectrum.
## Read in wavelengths
wave=as.numeric(sub(",","",read.table(x)[,1]))[1:3648]
wavekeep=wave>=range[1]&wave<=range[2]
wave=wave[wavekeep]
if(justwavelength) return(wave) #Return just wavelengths if that's all that's needed
## Read in Sample Values. For some reason the .jdx file says "% Reflection" but actually records
## just the sample values, so they have to be converted to % Reflection
sample=as.numeric(read.table(x)[,2][1:3648])[wavekeep]
## Read in the 'dark' reference values
dark=as.numeric(read.table(x,skip=28+3648+24)[,2][1:3648])[wavekeep]
## Test for a "bad dark" indicated by large SD
if(sd(dark)>50&darkstandard=="")
stop(paste("file: ",x,"had a dark measurement with an SD>50, indicating something is wrong. Please specify an alternative"))
if(sd(dark)>50&length(darkstandard)>1) {
dark=as.numeric(read.table(darkstandard,skip=28+3648+24)[,2][1:3648])[wavekeep]
print(paste("file: ",f,"had a bad dark, using darkstandard specified in function call, using alternative file:",darkstandard))}
## Read in the white reference values
white=as.numeric(read.table(x,skip=28+3648+24+3648+24)[,2][1:3648])[wavekeep]
## If an alternative white reference was indicated, use it. This is useful if you know the white reference is
## faulty and you want to use one from another sample
if(length(whitestandard)>1) white=as.numeric(read.table(whitestandard,skip=28+3648+24+3648+24)[,2][1:3648])[wavekeep]
## Check data type from .jdx file to confirm that it is Reflectance or Intensity
type=as.character(scan(x,what="character",skip=18,nlines=1,quiet=T))[2]
if(!type%in%c("Reflection","Intensity")) stop("File not Reflectance or Intensity, check data")
## Calculate % Reflectance
reflection=((sample-dark)/(white-dark))*100
## Set bounds on reflectance values. Often in the low sensitivity regions (<450 and >950, there are some strange values)
if(clip){
reflection=ifelse(reflection>0,reflection,0)
reflection=ifelse(reflection<100,reflection,100)
}
## build the output
filename=strsplit(x,"/")[[1]]
spectra=data.frame(temp=reflection)
colnames(spectra)=sub(".jdx","",filename[length(filename)])
rownames(spectra)=wave
print(paste("Finished importing file: ",x))
return(spectra)
}
view raw gistfile1.r hosted with ❤ by GitHub


You can use this to process a folder of spectra and return a data frame with a column for each file as follows:

sf=list.files("Spectra",pattern="*.jdx",full=T)
s=do.call(cbind,lapply(sf, spectracompile))
view raw gistfile1.r hosted with ❤ by GitHub


And here’s another function that will calculate NDVI similar to several satellite products:
calcndvi=function(x,wave=wave,sensor="MODIS") {
## calculate NDVI from reflectance vector
if(sensor=="AVHRR") {red=c(550,680); ired=c(730,1000)}
if(sensor=="TM") {red=c(630,690); ired=c(760,900)}
if(sensor=="MODIS") {red=c(620,670); ired=c(841,876)}
redr=which(wave>=red[1]&wave<=red[2])
iredr=which(wave>=ired[1]&wave<=ired[2])
return((mean(x[iredr])-mean(x[redr]))/(mean(x[iredr])+mean(x[redr])))
}
view raw gistfile1.r hosted with ❤ by GitHub


You can apply this to the output of spectracompile() like this:
ndvi=apply(s,2,calcndvi,wave=as.numeric(rownames(s)),sensor="MODIS")
view raw gistfile1.r hosted with ❤ by GitHub



Here’s a few random spectra from various fynbos species:

To leave a comment for the author, please follow the link and comment on their blog: PlanetFlux.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)