HyperocR: Hyperspectral OCR data processing

Documented in read.hocr.L2.SAS

#'
#'@title Reads HyperOCR data collected using Satlantic
#'       Hyperspectral Surface Aquisition System (HyperSAS)
#'
#'@description
#'The function reads ASCII files of Level 2 (calibrated) produced by ProSoft
#'software and returns a large list contaning all data.
#'Internal darks measurements may be subtracted or not using the
#'logical parameter APPLY.DARK.
#'
#'@param fn is the file name
#'@param VERBOSE is a logical parameter indicating whether or
#'               not the reading steps are printed. Default is TRUE.
#'@param APPLY.DARK is a logical parameter indicating whether or not
#'                  the darks are subtracted to the measurements.
#'                  Default is TRUE.
#'@return It returns a large list with all the parameters extracted
#'from the L2 file.
#'
#'@details
#'The dark measurements recorded during the measurements are available
#'in the L2 file and could be apply to the measurements if requested.
#'The lat/lon position and time are extracted from the header file.
#'Sun elevation is calculated from those information and stored in the "dd"
#'data frame returned in the list.
#'
#'(NOTE: Only *.dat file generated by Prosoft version 7.7.16 are supported)
#
#'@seealso
#'
#'@author
#'      Simon Bélanger
#'
#'@examples
#'
#'
#'@export
#'@name read.hocr.L2.SAS
read.hocr.L2.SAS <- function(fn, VERBOSE=TRUE, APPLY.DARK=TRUE){

  if (VERBOSE) print(paste("Reading:", fn))

  # Opens the file
  id = file(fn, "r")

  # Reads the first header line
  line = unlist(strsplit(readLines(con=id, n =1), " "))

  # extract Prosoft version
  ProSoftVersion <- line[length(line)]
  if (VERBOSE) print(paste("ProsSoft Version:", ProSoftVersion))

  # counts the number of header lines
  nrec = 1
  while (line[1] != "character(0)"){
    line = strsplit(readLines(con=id, n =1), " ")
    nrec <- nrec+1
    #print(line[1])
    if (line != "character(0)") {
      if (line[[1]][1] == "CAL_FILE_NAMES") {
        calfiles <- unlist(strsplit(line[[1]][3],","))
        SATTHS   <- any(str_detect(calfiles, pattern = "SATTHS")) # Create a logical variable indicating if a THS sensor was present
        ncal = length(calfiles)
        print(paste("Number of calibration files is: ", ncal))
        print(calfiles)
        # Extract the first three charcters of the cal name
        Cal.Ids =  toupper(str_sub(calfiles,1,3))
        Block.Type = str_sub(Cal.Ids, 3,3)
        id.Dark = which(Block.Type =="D")
        id.Meas = which(Block.Type !="D")
      }
    }

  }
  nHeaderLines = nrec + 1
  if (VERBOSE) print(paste("Number of header lines:", nHeaderLines))

  # count the number of dark lines for Ed (or ES)
  nrec = 0
  line = strsplit(readLines(con=id, n =1), " ")
  while (line != "character(0)"){
    line = strsplit(readLines(con=id, n =1), " ")
    nrec <- nrec+1
  }
  nDarkLinesEd = nrec + 1
  if (VERBOSE) print(paste("Number of dark records for Ed:",nDarkLinesEd))

  # count the number of dark lines for Li
  nrec = 0
  line = strsplit(readLines(con=id, n =1), " ")
  while (line != "character(0)"){
    line = strsplit(readLines(con=id, n =1), " ")
    nrec <- nrec+1
  }
  nDarkLinesLi = nrec + 1
  if (VERBOSE) print(paste("Number of dark records for Li:",nDarkLinesLi))

  # count the number of dark lines for Li
  nrec = 0
  line = strsplit(readLines(con=id, n =1), " ")
  while (line != "character(0)"){
    line = strsplit(readLines(con=id, n =1), " ")
    nrec <- nrec+1
  }
  nDarkLinesLt = nrec + 1
  if (VERBOSE) print(paste("Number of dark records for Lt:",nDarkLinesLt))

  # count the number of Ed measurements
  line = strsplit(readLines(con=id, n =1), " ")
  nrec = 0
  while (line != "character(0)"){
    line = strsplit(readLines(con=id, n =1), " ")
    nrec <- nrec+1
  }
  nEdLines = nrec
  if (VERBOSE) print(paste("Number of Ed records:",nEdLines))

  # count the number of Li measurements
  line = strsplit(readLines(con=id, n =1), " ")
  nrec = 0
  while (line != "character(0)"){
    line = strsplit(readLines(con=id, n =1), " ")
    nrec <- nrec+1
  }
  nLiLines = nrec
  if (VERBOSE) print(paste("Number of Li records:",nLiLines))

  # count the number of Lt measurements
  line = strsplit(readLines(con=id, n =1), " ")
  nrec = 0
  while (line != "character(0)"){
    line = strsplit(readLines(con=id, n =1), " ")
    nrec <- nrec+1
  }
  nLtLines = nrec
  if (VERBOSE) print(paste("Number of Lt records:",nLtLines))

  if (SATTHS) {
    # count the number of Ancillary measurements
    line = strsplit(readLines(con=id, n =1), " ")
    nrec = 0
    while (line != "character(0)"){
      line = strsplit(readLines(con=id, n =1), " ")
      nrec <- nrec+1
    }
    nAncLines = nrec
    if (VERBOSE) print(paste("Number of Ancillary records:",nAncLines))
    close(id)
  } else {
    if (VERBOSE) print("No Tilt Heading Sensor (THS) detected")
    close(id)
  }


  # Reads and store the header
  id = file(fn, "r")
  Header = rep("NA", nHeaderLines-1)
  for (i in 1:(nHeaderLines-1)) {
    Header[i] = readLines(con=id, n =1)
  }
  close(id)

  if (ProSoftVersion == "7.7.16_6" | ProSoftVersion == "7.7.19_2") {
    # Reads Ed darks and extract the wavelenght and time
    df.Ed.Darks =  read.table(fn, skip = nHeaderLines, nrows = (nDarkLinesEd-3), header=T)
    XLambda = names(df.Ed.Darks)[3:139]
    Ed.wl = as.numeric(str_sub(XLambda, 2,7))

    Ed.Darks = as.matrix(df.Ed.Darks[,3:139])
    DateDay = df.Ed.Darks[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    Hour = df.Ed.Darks[,147]

    Ed.Darks.Time = as.POSIXct(paste(Year,DOY,Hour),
                               format="%Y %j %H:%M:%S",tz="GMT")

    # Reads Li darks and extract the wavelenght and time
    df.Li.Darks =  read.table(fn, skip = nHeaderLines+nDarkLinesEd, nrows = (nDarkLinesLi-3), header=T)
    XLambda = names(df.Li.Darks)[3:139]
    Li.wl = as.numeric(str_sub(XLambda, 2,7))

    Li.Darks = as.matrix(df.Li.Darks[,3:139])
    if (ProSoftVersion == "7.7.19_2") DateDay = df.Li.Darks[,145] else DateDay = df.Li.Darks[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    if (ProSoftVersion == "7.7.19_2") Hour = df.Li.Darks[,146] else Hour = df.Li.Darks[,147]

    Li.Darks.Time = as.POSIXct(paste(Year,DOY,Hour),
                               format="%Y %j %H:%M:%S",tz="GMT")

    # Reads Lt darks and extract the wavelenght and time
    df.Lt.Darks =  read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi,
                              nrows = (nDarkLinesLt-3), header=T)
    XLambda = names(df.Lt.Darks)[3:139]
    Lt.wl = as.numeric(str_sub(XLambda, 2,7))

    Lt.Darks = as.matrix(df.Lt.Darks[,3:139])
    if (ProSoftVersion == "7.7.19_2") DateDay = df.Lt.Darks[,145] else DateDay = df.Lt.Darks[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    if (ProSoftVersion == "7.7.19_2") Hour = df.Lt.Darks[,146] else Hour = df.Lt.Darks[,147]

    Lt.Darks.Time = as.POSIXct(paste(Year,DOY,Hour),
                               format="%Y %j %H:%M:%S",tz="GMT")


    # Read the Ed data
    df.Ed = read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+nDarkLinesLt,
                       nrows = (nEdLines-3), header=T)

    Ed = as.matrix(df.Ed[,3:139])
    DateDay = df.Ed[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    Hour = df.Ed[,147]
    Ed.Time = as.POSIXct(paste(Year,DOY,Hour),
                         format="%Y %j %H:%M:%S",tz="GMT")
    if (VERBOSE) print("Ed")
    # Read the Li data
    df.Li = read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+nDarkLinesLt+nEdLines+1,
                       nrows = (nLiLines-3), header=T)

    Li = as.matrix(df.Li[,3:139])
    if (ProSoftVersion == "7.7.19_2") DateDay = df.Li[,145] else DateDay = df.Li[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    if (ProSoftVersion == "7.7.19_2") Hour = df.Li[,146] else Hour = df.Li[,147]
    Li.Time = as.POSIXct(paste(Year,DOY,Hour),
                         format="%Y %j %H:%M:%S",tz="GMT")
    if (VERBOSE) print("Li")

    # Read the Lt data
    df.Lt = read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+nDarkLinesLt+nEdLines+nLiLines+2,
                       nrows = (nLtLines-2), header=T)

    Lt = as.matrix(df.Lt[,3:139])
    if (ProSoftVersion == "7.7.19_2") DateDay = df.Lt[,145] else DateDay = df.Lt[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    if (ProSoftVersion == "7.7.19_2") Hour = df.Lt[,146] else Hour = df.Lt[,147]
    Lt.Time = as.POSIXct(paste(Year,DOY,Hour),
                         format="%Y %j %H:%M:%S",tz="GMT")
    if (VERBOSE) print("Lt")

    # compute mean darks
    mean.Ed.Dark = apply(Ed.Darks,2,mean, na.rm=T)
    mean.Lt.Dark = apply(Lt.Darks,2,mean, na.rm=T)
    mean.Li.Dark = apply(Li.Darks,2,mean, na.rm=T)

    # Apply DARK to measurements
    if (APPLY.DARK) {
      if (VERBOSE) print("Apply internal darks to Ed, Li and Lt measurements")

      mean.Ed.Dark.m = matrix(mean.Ed.Dark, ncol=length(mean.Ed.Dark), nrow=length(Ed[,1]), byrow = T)
      mean.Lt.Dark.m = matrix(mean.Lt.Dark, ncol=length(mean.Lt.Dark), nrow=length(Lt[,1]), byrow = T)
      mean.Li.Dark.m = matrix(mean.Li.Dark, ncol=length(mean.Li.Dark), nrow=length(Li[,1]), byrow = T)

      Ed = Ed - mean.Ed.Dark.m
      Lt = Lt - mean.Lt.Dark.m
      Li = Li - mean.Li.Dark.m


    } else {
      if (VERBOSE) print("WARNING: DARKS ARE NOT APPLY")

    }

    # Read Ancillary
    if (SATTHS) {
      df.Anc = read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+nDarkLinesLt
                          +nEdLines+nLiLines+nLtLines+3,
                          nrows = (nAncLines), header=T)

      names(df.Anc) = c("FRAME" , "TIMER" ,"ROLL", "PITCH",     "MAG_X", "MAG_Y",  "MAG_Z" ,
                        "COMP" , "DATETAG",   "TIME" ,"COUNTS")

      DateDay = df.Anc[,9]
      Year = as.numeric(str_sub(DateDay, 1,4))
      DOY = as.numeric(str_sub(DateDay, 5,7))
      Hour = df.Anc[,10]
      Anc.Time = as.POSIXct(paste(Year,DOY,Hour),
                            format="%Y %j %H:%M:%S",tz="GMT")
    } else {
      print("No THS to read")
    }




  }

  if (ProSoftVersion == "7.7.9") {

    # Reads and extracts Ed wavelengths
    Ed.head <- read.table(fn, skip = nHeaderLines, nrows = 1, header=F,comment.char = "!")
    Xlambda <- data.frame(Ed.head[which(as.character(Ed.head)!=1)],row.names=NULL)
    Ed.wl <- as.numeric(Xlambda)

    # Reads Ed darks and time
    df.Ed.Darks =  read.table(fn, skip = nHeaderLines,
                              nrows = (nDarkLinesEd-3), header=F)
    Ed.Darks = as.matrix(df.Ed.Darks[,3:139])
    DateDay = df.Ed.Darks[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    Hour = df.Ed.Darks[,147]

    Ed.Darks.Time = as.POSIXct(paste(Year,DOY,Hour),
                               format="%Y %j %H:%M:%S",tz="GMT")

    # Reads and extracts Li wavelengths
    Li.head <- read.table(fn, skip = nHeaderLines+nDarkLinesEd, nrows = 1, header=F)
    Xlambda <- data.frame(Li.head[which(as.character(Li.head)!=1)],row.names=NULL)
    Li.wl = as.numeric(Xlambda)

    # Reads Li darks and extract the wavelenght and time
    df.Li.Darks =  read.table(fn, skip = nHeaderLines+nDarkLinesEd+1,
                              nrows = (nDarkLinesLi-3), header=F)
    Li.Darks = as.matrix(df.Li.Darks[,3:139])
    DateDay = df.Li.Darks[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    Hour = df.Li.Darks[,147]

    Li.Darks.Time = as.POSIXct(paste(Year,DOY,Hour),
                               format="%Y %j %H:%M:%S",tz="GMT")

    # Reads and extracts Lt wavelengths
    Lt.head <- read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi,
                          nrows = 1, header=F)
    Xlambda <- data.frame(Lt.head[which(as.character(Lt.head)!=1)],row.names=NULL)
    Lt.wl = as.numeric(Xlambda)

    # Reads Lt darks and extract the wavelenght and time
    df.Lt.Darks =  read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+1,
                              nrows = (nDarkLinesLt-3), header=F)
    Lt.Darks = as.matrix(df.Lt.Darks[,3:139])
    DateDay = df.Lt.Darks[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    Hour = df.Lt.Darks[,147]

    Lt.Darks.Time = as.POSIXct(paste(Year,DOY,Hour),
                               format="%Y %j %H:%M:%S",tz="GMT")

    # Read the Ed data
    df.Ed = read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+nDarkLinesLt+1,
                       nrows = (nEdLines-2), header=F)
    Ed = as.matrix(df.Ed[,3:139])
    DateDay = df.Ed[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    Hour = df.Ed[,147]
    Ed.Time = as.POSIXct(paste(Year,DOY,Hour),
                         format="%Y %j %H:%M:%S",tz="GMT")
    if (VERBOSE) print("Ed")
    # Read the Li data
    df.Li = read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+nDarkLinesLt+nEdLines+2,
                       nrows = (nLiLines-2), header=F)

    Li = as.matrix(df.Li[,3:139])
    DateDay = df.Li[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    Hour = df.Li[,147]
    Li.Time = as.POSIXct(paste(Year,DOY,Hour),
                         format="%Y %j %H:%M:%S",tz="GMT")
    if (VERBOSE) print("Li")

    # Read the Lt data
    df.Lt = read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+nDarkLinesLt+nEdLines+nLiLines+3,
                       nrows = (nLtLines-2), header=F)

    Lt = as.matrix(df.Lt[,3:139])
    DateDay = df.Lt[,146]
    Year = as.numeric(str_sub(DateDay, 1,4))
    DOY = as.numeric(str_sub(DateDay, 5,7))
    Hour = df.Lt[,147]
    Lt.Time = as.POSIXct(paste(Year,DOY,Hour),
                         format="%Y %j %H:%M:%S",tz="GMT")
    if (VERBOSE) print("Lt")

    # compute mean darks
    mean.Ed.Dark = apply(Ed.Darks,2,mean, na.rm=T)
    mean.Lt.Dark = apply(Lt.Darks,2,mean, na.rm=T)
    mean.Li.Dark = apply(Li.Darks,2,mean, na.rm=T)

    # Apply DARK to measurements
    if (APPLY.DARK) {
      if (VERBOSE) print("Apply internal darks to Ed, Li and Lt measurements")

      mean.Ed.Dark.m = matrix(mean.Ed.Dark, ncol=length(mean.Ed.Dark), nrow=length(Ed[,1]), byrow = T)
      mean.Lt.Dark.m = matrix(mean.Lt.Dark, ncol=length(mean.Lt.Dark), nrow=length(Lt[,1]), byrow = T)
      mean.Li.Dark.m = matrix(mean.Li.Dark, ncol=length(mean.Li.Dark), nrow=length(Li[,1]), byrow = T)

      Ed = Ed - mean.Ed.Dark.m
      Lt = Lt - mean.Lt.Dark.m
      Li = Li - mean.Li.Dark.m


    } else {
      if (VERBOSE) print("WARNING: DARKS ARE NOT APPLY")

    }

    if (SATTHS) {
      # Read Ancillary
      df.Anc = read.table(fn, skip = nHeaderLines+nDarkLinesEd+nDarkLinesLi+nDarkLinesLt
                          +nEdLines+nLiLines+nLtLines+4,
                          nrows = (nAncLines), header=F)

      names(df.Anc) = c("FRAME" , "TIMER" ,"ROLL", "PITCH",     "MAG",
                        "COMP" , "DATETAG",   "TIME" ,"COUNTS")

      DateDay = df.Anc[,7]
      Year = as.numeric(str_sub(DateDay, 1,4))
      DOY = as.numeric(str_sub(DateDay, 5,7))
      Hour = df.Anc[,8]
      Anc.Time = as.POSIXct(paste(Year,DOY,Hour),
                            format="%Y %j %H:%M:%S",tz="GMT")
    } else {
      if (VERBOSE) print("No THS to read")
    }




  }

  # Extract information from header
  nHeaderLines = length(Header) - 1
  for (i in 1:nHeaderLines) {
    if (str_sub(Header[i],1,8) == "LATITUDE"){
      if (VERBOSE) print("Extract Latitude")
      split = str_split(Header[i], " ")
      deg = as.numeric(unlist(split)[3])
      if (str_length(unlist(split)[4]) == 6) {
        minutes = as.numeric(str_sub(unlist(split)[4],1,5))
      }
      if (str_length(unlist(split)[4]) == 7) {
        minutes = as.numeric(str_sub(unlist(split)[4],1,6))
      }
      lat = deg + (minutes/60)
      if (unlist(split)[5] == "N") {lat=lat}else{lat=-lat}
      if (VERBOSE) print(lat)
    }

    if (str_sub(Header[i],1,9) == "LONGITUDE"){
      if (VERBOSE) print("Extract Longitude")
      split = str_split(Header[i], " ")
      deg = as.numeric(unlist(split)[3])
      if (str_length(unlist(split)[4]) == 6) {
        minutes = as.numeric(str_sub(unlist(split)[4],1,5))
      }
      if (str_length(unlist(split)[4]) == 7) {
        minutes = as.numeric(str_sub(unlist(split)[4],1,6))
      }

      lon = deg + (minutes/60)
      if (unlist(split)[5] == "E") {lon=lon}else{lon=-lon}
      if (VERBOSE) print(lon)
    }

    if (str_sub(Header[i],1,4) == "ZONE") {
      if (VERBOSE) print("Extract Time Zone")

      if (str_sub(Header[i],8,10) == "Mis") {
        if (VERBOSE) print("Time zone is missing, GMT set as default")
        tz = "GMT"
        } else {
        tz = str_sub(Header[i],8,10)
        if (VERBOSE) print(tz)
        }
    }

    if (str_sub(Header[i],1,10) == "TIME-STAMP"){
      if (VERBOSE) print("Extract Time Stamp")
      Sys.setlocale("LC_TIME", "en_CA.UTF-8")
      TimeStamp = str_sub(Header[i],18,42)
      dte = as.POSIXct(as.character(TimeStamp), format="%b %d %H:%M:%S %Y", tz=tz)
      if (VERBOSE) print(dte)
    }
  }

  # Derive data from above information
  # sun-zenith angle
  day <- day(dte)
  month <- month(dte)
  year <- year(dte)
  hour <- hour(dte)
  minute <- minute(dte)
  second <- second(dte)
  ah <- hour + minute / 60 + second / 3600

  sunpos <- possol(month, day, ah, lon, lat)
  sunzen <- sunpos[1]
  sunazim <-sunpos[2]

  dd = list(longitude=lon, latitude=lat, sunzen = sunzen, sunazim=sunazim, date=dte)


  if (SATTHS) {
    SAS = list(Header=Header, Ed.wl=Ed.wl, Li.wl = Li.wl, Lt.wl = Lt.wl,
               Ed=Ed, Li=Li, Lt=Lt,
               Ed.Darks=Ed.Darks, Li.Darks=Li.Darks, Lt.Darks=Lt.Darks,
               Ed.Dark.mean=mean.Ed.Dark, Li.Dark.mean=mean.Li.Dark, Lt.Dark.mean=mean.Lt.Dark,
               Ed.Time =Ed.Time, Li.Time =Li.Time, Lt.Time = Lt.Time, Anc.Time = Anc.Time,
               Ed.Darks.Time = Ed.Darks.Time, Li.Darks.Time = Li.Darks.Time, Lt.Darks.Time = Lt.Darks.Time,
               Anc = df.Anc,
               dd=dd)
  } else {
    SAS = list(Header=Header, Ed.wl=Ed.wl, Li.wl = Li.wl, Lt.wl = Lt.wl,
               Ed=Ed, Li=Li, Lt=Lt,
               Ed.Darks=Ed.Darks, Li.Darks=Li.Darks, Lt.Darks=Lt.Darks,
               Ed.Dark.mean=mean.Ed.Dark, Li.Dark.mean=mean.Li.Dark, Lt.Dark.mean=mean.Lt.Dark,
               Ed.Time =Ed.Time, Li.Time =Li.Time, Lt.Time = Lt.Time, Anc.Time = NA,
               Ed.Darks.Time = Ed.Darks.Time, Li.Darks.Time = Li.Darks.Time, Lt.Darks.Time = Lt.Darks.Time,
               Anc = NA,
               dd=dd)
  }


  return(SAS)

}