R script to correct seismic file time stamps and structure

## PART 0 - SCRIPT INFORMATION SECTION ----------------------------------------

#' Script to correct the file name time stamp (and file type) of corrupt files
#' 
#' This script was written to correct five minute long mseed file snippets that 
#' were recorded by a 6-channel Centaur data logger. The file snippets were 
#' copied from the Centaur by a ssh command, which was triggered by a crontab 
#' on a Raspberry Pi every 5 minutes. The job resulted in about one second time 
#' delay and sometimes too short msseed files. In addition, not every 5 minute
#' snippet was actually recorded due to severe weather and data telemetry 
#' constraints at the deployment site. Finally, the time stamp of the mseed 
#' files does not correspond to the data they contain but the time the file 
#' was created. Therefore it became necessary to fix the corrupted files.
#' 
#' This script reads all available mseed files via ObsPy, splits the six 
#' channels into separate signal vectors, parses the time stamps and station 
#' IDs for each signal vector, creates correct eseis objects, and saves the 
#' data as sac files in an appropriate directory structure. 
#' 
#' The script was written as a patch solver. It should be used with care and 
#' the scientist should be aware of the consequences of the data manipulation. 
#' 
#' Author: Michael Dietze, Section 4.6 GFZ Potsdam (mdietze@gfz-potsdam.de)
#' 
#' Script version: 0.1.0, 05 Spetember 2019
#' 
#' Changes to previous versions
#'   - not applicable
#' 
#' Requirements & dependencies: 
#'   - R 3.4.4
#'   - eseis 0.5.0
#'   - reticulate 1.13
#'   
#'   - ObsPy 

## PART 1 - DEFINITION OF OBJECTS AND IMPORT OF DATA --------------------------

## load package
library(eseis)

## make ObsPy available
o <- reticulate::import("obspy")

## load station info file
stations <- read.table(file = paste0("data/stations/station_info_01.txt"), 
                       header = TRUE, 
                       stringsAsFactors = FALSE)

## set mseed directory
dir_mseed = paste0("data/seismic/centaur/service_01/")

## set output directory
dir_output <- "data/seismic/sac_redate/"

## list available mseed files
files_mseed <- list.files(path = dir_mseed,
                          recursive = TRUE,
                          full.names = TRUE, 
                          pattern = ".mseed")

## PART 2 - PROCESSING OF ALL FILES IN A LOOP ---------------------------------

for(i in 1:length(files_mseed)) {
  
  ## read and prepare files
  s <- try(eseis::aux_obspyeseis(o$read(files_mseed[i])))
  s <- try(s[order(names(s))])
  
  ## process channel-wise
  try(for(j in 1:length(s)) {
    
    ## assign channel to work with
    s_j <- s[[j]]
    
    ## extract start time
    s_start <- try(s_j$meta$starttime)
    
    ## assign station ID
    ID_out <- try(stations$ID[match(x = names(s)[j], 
                                    table = stations$Centaur_chn)])
    
    ## create output file name
    name_out <- try(paste0(dir_output,
                           format(x = s_start, format = "%Y/%j/"),
                           ID_out,
                           ".",
                           format(x = s_start, format = "%y.%j.%H.%M.%S"),
                           ".BHZ.SAC"))
    
    ## check/create directory year
    try(if(dir.exists(paths = paste0(dir_output,
                                     format(x = s_start, 
                                            format = "%Y"))) == FALSE) {
      
      dir.create(paste0(dir_output,
                        format(x = s_start, format = "%Y")))
    })
    
    ## check/create directory JD
    try(if(dir.exists(paths = paste0(dir_output,
                                     format(x = s_start, 
                                            format = "%Y/%j"))) == FALSE) {
      
      dir.create(paste0(dir_output,
                        format(x = s_start, format = "%Y/%j")))
    })
    
    ## assign meta data
    s_j$signal <- as.numeric(s_j$signal)
    s_j$meta$sensor <- "PE6B"
    s_j$meta$logger <- "Centaur"
    s_j$meta$gain <- 40
    s_j$meta$station <- ID_out
    s_j$meta$latitude <- 50
    s_j$meta$longitude <- 13
    s_j$meta$elevation <- 1000
    s_j$meta$depth <- 0
    
    ## write output SAC file
    write_sac(data = s_j, 
              file = name_out, 
              unit = 1)
  })
  
  ## print progress
  print(paste(i, "of", length(files_mseed)))
}