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.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | |
} |
You can use this to process a folder of spectra and return a data frame with a column for each file as follows:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
sf=list.files("Spectra",pattern="*.jdx",full=T) | |
s=do.call(cbind,lapply(sf, spectracompile)) |
And here’s another function that will calculate NDVI similar to several satellite products:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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]))) | |
} |
You can apply this to the output of spectracompile() like this:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
ndvi=apply(s,2,calcndvi,wave=as.numeric(rownames(s)),sensor="MODIS") |
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.