R/echoIBM.calibrate.R
In arnejohannesholmin/echoIBM: Simulation of multibeam sonar data
Defines functions
echoIBM.calibrate
Documented in
echoIBM.calibrate
#*********************************************
#*********************************************
#' Calibrates (and writes the calibration data to file) the simulation of observations from the MS70 sonar performed by echoIBM(). Not used in any of the other echoIBM code files.
#'
#' @param directory the path to the directory in which to put the calibration file (use TRUE for the default directory file.path(echoIBM_datasets_directory(), "Resources", "Calibration")).
#' @param cruise is either the idenfication number of the cruise, given as specified by the IMR (yyyynnn), or the path to the directory containing the event on which the calibration should be based.
#' @param event is the identifier of the event, either given as the number of the event, a string contained in the name of the event, or the path of the event directory.
#' @param esnm is the name of the acoustical instrument, given as a four character string. See sonR_implemented() for the implemented systems. May be given in 'data', and in lower case.
#' @param type is "sv" for calibrating the simulation model with respect to volume scattering 'sv' and "ts" for calibrating the simulation model with respect to backscattering cross section / target strength.
#' @param runs is an integer giving the number of simulated pings of which the mean is taken as the final calibration factors.
#' @param add is the angle in degrees to add on either end of the insonifying volume of the acoustical instrument, defining the volume in which point targets are distributed.
#' @param N is the number of point targets used in the calibration.
#' @param j is the sampling interval to use for the calibration (should be larger than 2).
#' @param sgbs is the backscattering cross section of the point targets.
#' @param max.memory is used in echoIBM to set the maximum memory for echoIBM to occupy.
#' @param filetag is a string to add to the file name.
#'
#' @return
#'
#' @examples
#' \dontrun{}
#'
#' @importFrom SimradRaw soundbeam_range
#' @importFrom sonR get.specs.esnm read.event runif.sph
#' @importFrom TSD NAs sph2car strff write.TSD zeros
#'
#' @export
#' @rdname echoIBM.calibrate
#'
echoIBM.calibrate <- function(directory=NULL, cruise="S2009116", event=1, esnm="MS70", type=c("sv","ts"), dire=NULL, runs=1, add=30, N=1e5, j=100, sgbs=1e-6, max.memory=1e9, method=c("closest","linear"), max.radius=0.2, seed=0, cores=1, usemean=FALSE, nameadd=""){
############ AUTHOR(S): ############
# Arne Johannes Holmin
############ LANGUAGE: #############
# English
############### LOG: ###############
# Start: 2010-03-20 - Clean version.
# Update: 2011-06-25 - Major changes: The option 'type' added, for choosing between calibrating for sv or for TS. A factor 1/(4*pi) added to the true 'sv'. File written to the directory 'directory', while the previous version only had the option to write to a specific file. The file name is as follows: paste("CalibrationTable_",toupper(esnm),"_Cruise_",cruise,"_sigma_bs.beams",sep="").
# Update: 2011-06-26 - Added the option 'runs', for running multiple simulations of fewer calibration targets (saving memory).
# Update: 2011-08-24 - Changed from using 'epsl', 'acca' and 'ssil' to using 'epss' and 'sgbs'. Tested using the script "Test echoIBM.calibrate.R".
# Update: 2012-02-13 - Changed the method radically, by using echoIBM exactly in the way an ordinary simulation would be performed, in order to keep the same control over changes in the simulation model without making the same changes in this function (the previous version used echoIBM.oneping.oneschool). The data are then read back in to the function, and the calibration factors calculated.
# Update: 2013-02-25 - Saved some minor changes several months ago.
# Last: 2015-02-22 - Added support for multiple tilt angles through 'dire'.
##################################################
##################################################
########## Preparation ##########
if(isTRUE(directory)){
directory <- file.path(echoIBM_datasets_directory(), "Resources", "Calibration")
}
# 'beams' must be given, either as a list of beam configuration variables, or as the path to a beams configuration file:
data <- read.event(var="beams",cruise=cruise,event=event,esnm=esnm)
# Warning if 'data' is empty:
if(length(data)<2){
stop(paste0("No beams data in cruise ", cruise, ", event ", event))
}
data$rres <- soundbeam_range(data, "rres")
if(j-2>min(data$lenb)){
stop(paste("The sampling interval index 'j' chosen too large for the minimum beam length ",min(data$lenb),sep=""))
}
# Get name of the acoustical instrument:
if(is.null(esnm)){
esnm <- data$esnm
}
esnm <- toupper(esnm)
# Get the number of beams:
Ni <- length(data$dira)
########## Execution and output ##########
# The calibration procedure of echoIBM for volume backscattering coefficient 'sv' is as follows: Draw 'N' uniformly distributed spherical targets of backscattering cross section 'sgbs' covering the volume affecting all voxels corresponding to sampling interval 'j' (see appendix 2 of the documentation of echoIBM). Define the true 'sv' to be N*sgbs/V, where 'V' is the size of the volume covered by sampling interval 'j' in the system of voxels. Then run the simulations and compare to the true 'sv' by defining calibration factors equal to true 'sv' / colSums(simulated 'sv').
if(strff("s",type[1])){
### (1) Create the directory of the calibration simulation: ###
#if(length(directory)>0 && !identical(directory,FALSE)){
calname <- paste("Calibration",esnm,"Cruise",cruise,"sv",sep="_")
#}
#else{
# calname <- paste("Calibration",esnm,"Cruise",cruise,"sv_TEST",sep="_")
#}
caldir <- file.path(tempdir(), "Calibration", "Events", calname, esnm, "tsd")
#caldir <- paste("~/Data/echoIBM/Calibration/Events/",calname,"/",esnm,"/tsd",sep="")
### if(isTRUE(file.info(caldir)$isdir)){
### ans=readline(paste("The calibration case (",caldir,") alreaddy exists. Replace? (y/n)"))
### if(!identical(ans,"y")){
### stop("Calibration terminated")
### }
### else{
### unlink(caldir, recursive = TRUE)
### warning("The existing calibration case replaced")
### }
### }
unlink(caldir, recursive=TRUE, force=TRUE)
dir.create(caldir,recursive=TRUE)
#if(is.na(file.info(caldir)$isdir)){
# dir.create(caldir,recursive=TRUE)
#}
#else{
# ans=readline(paste("The calibration case (",caldir,") alreaddy exists. Replace? (y/n)"))
# if(!identical(ans,"y")){
# stop("Calibration terminated")
# }
# else{
# tempdir=paste("~/Data/echoIBM/Calibration/Events/",calname,"/",esnm,"/temp",sep="")
# if(isTRUE(file.info(tempdir)$isdir)){
# unlink(tempdir, recursive = TRUE)
# }
# warning("The existing calibration case replaced")
# }
#}
### (2) Write the beam-configuration data: ###
cat("Write the beam-configuration data...\n")
write.TSD(data, paste(caldir,"/",calname,".beams",sep=""), numt=1)
### (3) Write the vessel data: ###
cat("Write the vessel data...\n")
vessel <- list(psxv=zeros(runs),psyv=zeros(runs),pszv=zeros(runs),rtzv=zeros(runs),utim=sequence(runs)-1)
write.TSD(vessel, paste(caldir,"/",calname,".vessel",sep=""), numt=runs)
### (4) Write the static school data: ###
schoolstatic <- list(pbpf="ps",sgbs=sgbs,obln=1,gamw=0,gaml=0)
cat("Write the static school data...\n")
schoolstaticFile <- paste(caldir,"/",calname,"_pointTargets_static.school",sep="")
write.TSD(schoolstatic, schoolstaticFile, numt=1)
# If 'dire' is given, run through its elements and calibrate for each value, equal for all beams:
if(length(dire)>0){
cali <- NAs(Ni, length(dire))
for(k in seq_along(dire)){
cat("Calibration for tilt angle ",k," (of",length(dire),") ...\n")
# Use the input elevation angle:
thisdire <- list(dire=rep(dire[k],length(data$lenb)))
### (5) Generate the target positions and write the school file: ###
region <- get.specs.esnm(c(thisdire, data), esnm=esnm, add=add)$region
#dr <- data$asps*data$sint/2
r0 <- (j-2)*data$rres
r1 <- j*data$rres
# Total volume in which the target positions are drawn. Also see appendix 2 of the documentation of echoIBM:
V <- vol.sph(r=c(r0,r1),theta=region[,1],phi=region[,2])$volx
# The true volume backscatter, equal to N*sigma_bs/V = N*A_0/V for spherical targets:
truesv <- sgbs*N/V
# Generate the target positions and write to the school file at each run:
cat("Write the dynamic school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_dynamic.school",sep="")
for(i in seq_len(runs)){
cat("Run ",i,":\t",sep="")
### Generate positions corresponding to uniformly distributed objects in space, from which to simulate observations from the MS70 sonar. See appendix 2 of the documentation of echoIBM: ###
xyz <- runif.sph(N,r=c(r0,r1),theta=region[,1],phi=region[,2],sph.out=FALSE)
#plot3d(xyz,size=0.3)
# Add transducer position and cut the positions above the sea level:
xyz[,3] <- xyz[,3]+data$psze[1]
xyz <- xyz[xyz[,3]<0,]
# Update 'N':
thisN <- nrow(xyz)
cat(thisN,"point sources drawn\n")
# Add the target information (all rotations left out because of points targets):
schooldynamic <- list(psxf=xyz[,1],psyf=xyz[,2],pszf=xyz[,3],utim=i-1)
# Write the dynamic school data to file:
if(i==1){
write.TSD(schooldynamic, schooldynamicFile, reserve=runs-1, numt=1)
}
else{
write.TSD(schooldynamic, schooldynamicFile, numt=1, append=TRUE)
}
}
### (6) Run the simulation: ###
echoIBM(caldir, adds=thisdire, t="all", TVG.exp=2, calibrate=FALSE, noise="", max.memory=max.memory, max.radius=max.radius, method=method, esnm=esnm, parlist=list(seed=seed), keep.temp=FALSE)
### (7) Read the simulated data: ###
sv <- zeros(max(data$lenb),length(data$lenb))
for(i in seq_len(runs)){
thissv <- read.event(event=caldir,t=i,var="vbsc")$vbsc
sv <- sv+thissv
}
### (8) Store the calibration factors of the current tilt: ###
cali[,k] <- truesv/sv[j,]*runs
}
### (9) Write and return the calibration data: ###
# Apply the mean of all beams if usemean==TRUE (omnidirectional sonars such as Simrad SX90):
if(usemean){
#out=c(list(cali=apply(cali,2,function(xx) rep(mean(xx),length(xx))), cal0=cali, grde=dire, ctim=unclass(Sys.time())), data[setdiff(names(data),c("dire","ftim","mtim","utim"))])
out <- c(list(cali=apply(cali,2,function(xx) rep(mean(xx),length(xx))), cal0=cali, grde=dire, ctim=unclass(Sys.time())))
}
else{
#out=c(list(cali=cali, grde=dire, ctim=unclass(Sys.time())), data[setdiff(names(data),c("dire","ftim","mtim","utim"))])
out <- c(list(cali=cali, grde=dire, ctim=unclass(Sys.time())))
}
# Write the calibration data:
if(length(directory)>0 && !identical(directory,FALSE)){
con <- file.path(directory, paste("CalibrationTable_", toupper(esnm), "_Cruise_", cruise, if(nchar(nameadd)>0) "_", nameadd, "_sv.beams", sep=""))
write.TSD(out, con=con, numt=1)
}
# Return the calibration data:
c(out, list(truesv=truesv, caldir=caldir))
}
else{
### (5) Generate the target positions and write the school file: ###
region <- get.specs.esnm(data, esnm=esnm, add=add)$region
#dr=data$asps*data$sint/2
r0 <- (j-2)*data$rres
r1 <- j*data$rres
# Total volume in which the target positions are drawn. Also see appendix 2 of the documentation of echoIBM:
V <- vol.sph(r=c(r0,r1),theta=region[,1],phi=region[,2])$volx
# The true volume backscatter, equal to N*sigma_bs/V = N*A_0/V for spherical targets:
truesv <- sgbs*N/V
# Generate the target positions and write to the school file at each run:
cat("Write the dynamic school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_dynamic.school",sep="")
for(i in seq_len(runs)){
cat("Run ",i,":\t",sep="")
### Generate positions corresponding to uniformly distributed objects in space, from which to simulate observations from the MS70 sonar. See appendix 2 of the documentation of echoIBM: ###
xyz <- runif.sph(N,r=c(r0,r1),theta=region[,1],phi=region[,2],sph.out=FALSE)
#plot3d(xyz,size=0.3)
# Add transducer position and cut the positions above the sea level:
xyz[,3] <- xyz[,3]+data$psze[1]
xyz <- xyz[xyz[,3]<0,]
# Update 'N':
thisN <- nrow(xyz)
cat(thisN,"point sources drawn\n")
# Add the target information (all rotations left out because of points targets):
schooldynamic <- list(psxf=xyz[,1],psyf=xyz[,2],pszf=xyz[,3],utim=i-1)
# Write the dynamic school data to file:
if(i==1){
write.TSD(schooldynamic, schooldynamicFile, reserve=runs-1, numt=1)
}
else{
write.TSD(schooldynamic, schooldynamicFile, numt=1, append=TRUE)
}
}
### (6) Run the simulation: ###
echoIBM(caldir, t="all", TVG.exp=2, calibrate=FALSE, noise="", max.memory=max.memory, max.radius=max.radius, method=method, esnm=esnm, parlist=list(seed=seed), keep.temp=FALSE)
### (7) Read the simulated data: ###
sv <- zeros(max(data$lenb),length(data$lenb))
for(i in seq_len(runs)){
thissv <- read.event(event=caldir,t=i,var="vbsc")$vbsc
sv <- sv+thissv
}
### (8) Write the calibration data if required: ###
calfact <- truesv/sv[j,]*runs
### (9) Write and return the calibration data: ###
# Apply the mean of all beams if usemean==TRUE (omnidirectional sonars such as Simrad SX90):
if(usemean){
#out=c(list(cali=rep(mean(calfact),length(calfact)), cal0=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
out <- c(list(cali=rep(mean(calfact),length(calfact)), cal0=calfact, ctim=unclass(Sys.time())))
}
else{
#out=c(list(cali=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
out <- c(list(cali=calfact, ctim=unclass(Sys.time())))
}
# Write the calibration data:
if(length(directory)>0 && !identical(directory,FALSE)){
con <- file.path(directory,paste("CalibrationTable_",toupper(esnm),"_Cruise_",cruise,if(nchar(nameadd)>0) "_",nameadd,"_sv.beams",sep=""))
write.TSD(out, con=con, numt=1)
}
# Return the calibration data:
c(out, list(truesv=truesv, caldir=caldir))
}
}
# The calibration procedure for target strength 'TS' is as follows: Place one target on the acoustic axis of each beam, and simulate the TS using type="ts" in echoIBM(). This includes 40*log*r for the TVG function instead of 20*log*r as for the sv. Compare the simulated 'sigma_bs' to the true 'sigma_bs' by defining calibration factors equal to true 'sigma_bs' / colSums(simulated 'sigma_bs').
else if(strff("t",type[1])){
### (1) Create the directory of the calibration simulation: ###
if(length(directory)>0 && !identical(directory,FALSE)){
calname <- paste("Calibration",esnm,"Cruise",cruise,"sigma_bs",sep="_")
}
else{
calname <- paste("Calibration",esnm,"Cruise",cruise,"sigma_bs_TEST",sep="_")
}
caldir <- paste("~/Data/echoIBM/Calibration/Events/",calname,"/",esnm,"/tsd",sep="")
if(is.na(file.info(caldir)$isdir)){
dir.create(caldir,recursive=TRUE)
}
else{
ans <- readline("The calibration case alreaddy exists. Replace? (y/n)")
if(!identical(ans,"y")){
stop("Calibration terminated")
}
else{
tempdir <- paste("~/Data/echoIBM/Calibration/Events/",calname,"/",esnm,"/temp",sep="")
if(isTRUE(file.info(tempdir)$isdir)){
unlink(tempdir, recursive = TRUE)
}
warning("The existing calibration case replaced")
}
}
### (2) Write the beam-configuration data: ###
cat("Write the beam-configuration data...\n")
write.TSD(data, paste(caldir,"/",calname,".beams",sep=""), numt=1)
### (3) Write the vessel data: ###
cat("Write the vessel data...\n")
vessel <- list(psxv=zeros(Ni),psyv=zeros(Ni),pszv=zeros(Ni),rtzv=zeros(Ni),utim=sequence(Ni)-1)
write.TSD(vessel,paste(caldir,"/",calname,".vessel",sep=""),numt=Ni)
### (4) Write the static school data: ###
schoolstatic <- list(pbpf="ps",sgbs=sgbs,obln=1,gamw=0,gaml=0)
cat("Write the static school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_static.school",sep="")
write.TSD(schoolstatic, schooldynamicFile, numt=1)
# If 'dire' is given, run through its elements and calibrate for each value, equal for all beams:
if(length(dire)>0){
cali <- NAs(Ni, length(dire))
for(k in seq_along(dire)){
cat("Calibration for tilt angle ",k," (of",length(dire),") ...\n")
# Use the input elevation angle:
thisdire <- list(dire=rep(dire[k],length(data$lenb)))
### (5) Generate the target positions and write the school file: ###
rthph <- cbind((j-1)*data$rres, data$dira, thisdire$dire)
xyz <- sph2car(rthph)
xyz[,3] <- xyz[,3]+data$psze[1]
schooldynamic <- list(psxf=xyz[,1],psyf=xyz[,2],pszf=xyz[,3],utim=sequence(Ni)-1)
cat("Write the dynamic school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_dynamic.school",sep="")
write.TSD(schooldynamic, schooldynamicFile, numt=Ni)
### (6) Run the simulation: ###
echoIBM(caldir, adds=thisdire, t="all", TVG.exp=4, calibrate=FALSE, noise="", max.memory=max.memory, method=method, max.radius=max.radius, esnm=esnm, parlist=list(seed=seed), keep.temp=FALSE)
### (7) Read the simulated data: ###
thiscali <- zeros(Ni)
for(i in seq_len(Ni)){
thissv <- read.event(event=caldir,t=i,var="vbsc")$vbsc
thiscali[i] <- sum(thissv[,i])
}
### (8) Store the calibration factors of the current tilt: ###
cali[,k] <- sgbs/thiscali
}
### (9) Write and return the calibration data: ###
# Apply the mean of all beams if usemean==TRUE (omnidirectional sonars such as Simrad SX90):
if(usemean){
#out=c(list(cali=apply(cali,2,function(xx) rep(mean(xx),length(xx))), cal0=cali, grde=dire, ctim=unclass(Sys.time())), data[setdiff(names(data),c("dire","ftim","mtim","utim"))])
out <- c(list(cali=apply(cali,2,function(xx) rep(mean(xx),length(xx))), cal0=cali, grde=dire, ctim=unclass(Sys.time())))
}
else{
#out=c(list(cali=cali, grde=dire, ctim=unclass(Sys.time())), data[setdiff(names(data),c("dire","ftim","mtim","utim"))])
out <- c(list(cali=cali, grde=dire, ctim=unclass(Sys.time())))
}
# Write the calibration data:
if(length(directory)>0 && !identical(directory,FALSE)){
con <- file.path(directory,paste("CalibrationTable_",toupper(esnm),"_Cruise_",cruise,if(nchar(nameadd)>0) "_",nameadd,"_sigma_bs.beams",sep=""))
write.TSD(out, con=con, numt=1)
}
# Return the calibration data:
c(out, list(truesv=truesv, caldir=caldir))
}
else{
### (5) Generate the target positions and write the school file: ###
rthph <- cbind((j-1)*data$rres, data$dira, data$dire)
xyz <- sph2car(rthph)
xyz[,3] <- xyz[,3]+data$psze[1]
schooldynamic <- list(psxf=xyz[,1],psyf=xyz[,2],pszf=xyz[,3],utim=sequence(Ni)-1)
cat("Write the dynamic school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_dynamic.school",sep="")
write.TSD(schooldynamic, schooldynamicFile, numt=Ni)
### (6) Run the simulation: ###
echoIBM(caldir, t="all", TVG.exp=4, calibrate=FALSE, noise="", max.memory=max.memory, method=method, max.radius=max.radius, esnm=esnm, parlist=list(seed=seed), keep.temp=FALSE, cores=cores)
### (7) Read the simulated data: ###
cali <- zeros(Ni)
for(i in seq_len(Ni)){
thissv <- read.event(event=caldir,t=i,var="vbsc")$vbsc
cali[i] <- sum(thissv[,i])
}
### (8) Write the calibration data if required: ###
calfact <- sgbs/cali
# Apply the mean of all beams if usemean==TRUE (omnidirectional sonars such as Simrad SX90):
if(usemean){
#out=c(list(cali=rep(mean(calfact),length(calfact)), cal0=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
out <- c(list(cali=rep(mean(calfact),length(calfact)), cal0=calfact, ctim=unclass(Sys.time())))
}
else{
#out=c(list(cali=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
#out=c(list(cali=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
out <- c(list(cali=calfact, ctim=unclass(Sys.time())))
}
# Write the calibration data:
if(length(directory)>0 && !identical(directory,FALSE)){
con <- file.path(directory,paste("CalibrationTable_",toupper(esnm),"_Cruise_",cruise,if(nchar(nameadd)>0) "_",nameadd,"_sigma_bs.beams",sep=""))
write.TSD(out, con=con, numt=1)
}
# Return the calibration data:
c(out, list(truesv=truesv, caldir=caldir))
}
}
# Else return an error if the calibration type is not recognized:
else{
stop("Invalid calibration type (must be one of \"sv\" or \"ts\")")
}
##################################################
##################################################
}
arnejohannesholmin/echoIBM documentation built on April 14, 2024, 11:37 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions echoIBM.calibrate
Documented in echoIBM.calibrate
#*********************************************
#*********************************************
#' Calibrates (and writes the calibration data to file) the simulation of observations from the MS70 sonar performed by echoIBM(). Not used in any of the other echoIBM code files.
#'
#' @param directory the path to the directory in which to put the calibration file (use TRUE for the default directory file.path(echoIBM_datasets_directory(), "Resources", "Calibration")).
#' @param cruise is either the idenfication number of the cruise, given as specified by the IMR (yyyynnn), or the path to the directory containing the event on which the calibration should be based.
#' @param event is the identifier of the event, either given as the number of the event, a string contained in the name of the event, or the path of the event directory.
#' @param esnm is the name of the acoustical instrument, given as a four character string. See sonR_implemented() for the implemented systems. May be given in 'data', and in lower case.
#' @param type is "sv" for calibrating the simulation model with respect to volume scattering 'sv' and "ts" for calibrating the simulation model with respect to backscattering cross section / target strength.
#' @param runs is an integer giving the number of simulated pings of which the mean is taken as the final calibration factors.
#' @param add is the angle in degrees to add on either end of the insonifying volume of the acoustical instrument, defining the volume in which point targets are distributed.
#' @param N is the number of point targets used in the calibration.
#' @param j is the sampling interval to use for the calibration (should be larger than 2).
#' @param sgbs is the backscattering cross section of the point targets.
#' @param max.memory is used in echoIBM to set the maximum memory for echoIBM to occupy.
#' @param filetag is a string to add to the file name.
#'
#' @return
#'
#' @examples
#' \dontrun{}
#'
#' @importFrom SimradRaw soundbeam_range
#' @importFrom sonR get.specs.esnm read.event runif.sph
#' @importFrom TSD NAs sph2car strff write.TSD zeros
#'
#' @export
#' @rdname echoIBM.calibrate
#'
echoIBM.calibrate <- function(directory=NULL, cruise="S2009116", event=1, esnm="MS70", type=c("sv","ts"), dire=NULL, runs=1, add=30, N=1e5, j=100, sgbs=1e-6, max.memory=1e9, method=c("closest","linear"), max.radius=0.2, seed=0, cores=1, usemean=FALSE, nameadd=""){
############ AUTHOR(S): ############
# Arne Johannes Holmin
############ LANGUAGE: #############
# English
############### LOG: ###############
# Start: 2010-03-20 - Clean version.
# Update: 2011-06-25 - Major changes: The option 'type' added, for choosing between calibrating for sv or for TS. A factor 1/(4*pi) added to the true 'sv'. File written to the directory 'directory', while the previous version only had the option to write to a specific file. The file name is as follows: paste("CalibrationTable_",toupper(esnm),"_Cruise_",cruise,"_sigma_bs.beams",sep="").
# Update: 2011-06-26 - Added the option 'runs', for running multiple simulations of fewer calibration targets (saving memory).
# Update: 2011-08-24 - Changed from using 'epsl', 'acca' and 'ssil' to using 'epss' and 'sgbs'. Tested using the script "Test echoIBM.calibrate.R".
# Update: 2012-02-13 - Changed the method radically, by using echoIBM exactly in the way an ordinary simulation would be performed, in order to keep the same control over changes in the simulation model without making the same changes in this function (the previous version used echoIBM.oneping.oneschool). The data are then read back in to the function, and the calibration factors calculated.
# Update: 2013-02-25 - Saved some minor changes several months ago.
# Last: 2015-02-22 - Added support for multiple tilt angles through 'dire'.
##################################################
##################################################
########## Preparation ##########
if(isTRUE(directory)){
directory <- file.path(echoIBM_datasets_directory(), "Resources", "Calibration")
}
# 'beams' must be given, either as a list of beam configuration variables, or as the path to a beams configuration file:
data <- read.event(var="beams",cruise=cruise,event=event,esnm=esnm)
# Warning if 'data' is empty:
if(length(data)<2){
stop(paste0("No beams data in cruise ", cruise, ", event ", event))
}
data$rres <- soundbeam_range(data, "rres")
if(j-2>min(data$lenb)){
stop(paste("The sampling interval index 'j' chosen too large for the minimum beam length ",min(data$lenb),sep=""))
}
# Get name of the acoustical instrument:
if(is.null(esnm)){
esnm <- data$esnm
}
esnm <- toupper(esnm)
# Get the number of beams:
Ni <- length(data$dira)
########## Execution and output ##########
# The calibration procedure of echoIBM for volume backscattering coefficient 'sv' is as follows: Draw 'N' uniformly distributed spherical targets of backscattering cross section 'sgbs' covering the volume affecting all voxels corresponding to sampling interval 'j' (see appendix 2 of the documentation of echoIBM). Define the true 'sv' to be N*sgbs/V, where 'V' is the size of the volume covered by sampling interval 'j' in the system of voxels. Then run the simulations and compare to the true 'sv' by defining calibration factors equal to true 'sv' / colSums(simulated 'sv').
if(strff("s",type[1])){
### (1) Create the directory of the calibration simulation: ###
#if(length(directory)>0 && !identical(directory,FALSE)){
calname <- paste("Calibration",esnm,"Cruise",cruise,"sv",sep="_")
#}
#else{
# calname <- paste("Calibration",esnm,"Cruise",cruise,"sv_TEST",sep="_")
#}
caldir <- file.path(tempdir(), "Calibration", "Events", calname, esnm, "tsd")
#caldir <- paste("~/Data/echoIBM/Calibration/Events/",calname,"/",esnm,"/tsd",sep="")
### if(isTRUE(file.info(caldir)$isdir)){
### ans=readline(paste("The calibration case (",caldir,") alreaddy exists. Replace? (y/n)"))
### if(!identical(ans,"y")){
### stop("Calibration terminated")
### }
### else{
### unlink(caldir, recursive = TRUE)
### warning("The existing calibration case replaced")
### }
### }
unlink(caldir, recursive=TRUE, force=TRUE)
dir.create(caldir,recursive=TRUE)
#if(is.na(file.info(caldir)$isdir)){
# dir.create(caldir,recursive=TRUE)
#}
#else{
# ans=readline(paste("The calibration case (",caldir,") alreaddy exists. Replace? (y/n)"))
# if(!identical(ans,"y")){
# stop("Calibration terminated")
# }
# else{
# tempdir=paste("~/Data/echoIBM/Calibration/Events/",calname,"/",esnm,"/temp",sep="")
# if(isTRUE(file.info(tempdir)$isdir)){
# unlink(tempdir, recursive = TRUE)
# }
# warning("The existing calibration case replaced")
# }
#}
### (2) Write the beam-configuration data: ###
cat("Write the beam-configuration data...\n")
write.TSD(data, paste(caldir,"/",calname,".beams",sep=""), numt=1)
### (3) Write the vessel data: ###
cat("Write the vessel data...\n")
vessel <- list(psxv=zeros(runs),psyv=zeros(runs),pszv=zeros(runs),rtzv=zeros(runs),utim=sequence(runs)-1)
write.TSD(vessel, paste(caldir,"/",calname,".vessel",sep=""), numt=runs)
### (4) Write the static school data: ###
schoolstatic <- list(pbpf="ps",sgbs=sgbs,obln=1,gamw=0,gaml=0)
cat("Write the static school data...\n")
schoolstaticFile <- paste(caldir,"/",calname,"_pointTargets_static.school",sep="")
write.TSD(schoolstatic, schoolstaticFile, numt=1)
# If 'dire' is given, run through its elements and calibrate for each value, equal for all beams:
if(length(dire)>0){
cali <- NAs(Ni, length(dire))
for(k in seq_along(dire)){
cat("Calibration for tilt angle ",k," (of",length(dire),") ...\n")
# Use the input elevation angle:
thisdire <- list(dire=rep(dire[k],length(data$lenb)))
### (5) Generate the target positions and write the school file: ###
region <- get.specs.esnm(c(thisdire, data), esnm=esnm, add=add)$region
#dr <- data$asps*data$sint/2
r0 <- (j-2)*data$rres
r1 <- j*data$rres
# Total volume in which the target positions are drawn. Also see appendix 2 of the documentation of echoIBM:
V <- vol.sph(r=c(r0,r1),theta=region[,1],phi=region[,2])$volx
# The true volume backscatter, equal to N*sigma_bs/V = N*A_0/V for spherical targets:
truesv <- sgbs*N/V
# Generate the target positions and write to the school file at each run:
cat("Write the dynamic school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_dynamic.school",sep="")
for(i in seq_len(runs)){
cat("Run ",i,":\t",sep="")
### Generate positions corresponding to uniformly distributed objects in space, from which to simulate observations from the MS70 sonar. See appendix 2 of the documentation of echoIBM: ###
xyz <- runif.sph(N,r=c(r0,r1),theta=region[,1],phi=region[,2],sph.out=FALSE)
#plot3d(xyz,size=0.3)
# Add transducer position and cut the positions above the sea level:
xyz[,3] <- xyz[,3]+data$psze[1]
xyz <- xyz[xyz[,3]<0,]
# Update 'N':
thisN <- nrow(xyz)
cat(thisN,"point sources drawn\n")
# Add the target information (all rotations left out because of points targets):
schooldynamic <- list(psxf=xyz[,1],psyf=xyz[,2],pszf=xyz[,3],utim=i-1)
# Write the dynamic school data to file:
if(i==1){
write.TSD(schooldynamic, schooldynamicFile, reserve=runs-1, numt=1)
}
else{
write.TSD(schooldynamic, schooldynamicFile, numt=1, append=TRUE)
}
}
### (6) Run the simulation: ###
echoIBM(caldir, adds=thisdire, t="all", TVG.exp=2, calibrate=FALSE, noise="", max.memory=max.memory, max.radius=max.radius, method=method, esnm=esnm, parlist=list(seed=seed), keep.temp=FALSE)
### (7) Read the simulated data: ###
sv <- zeros(max(data$lenb),length(data$lenb))
for(i in seq_len(runs)){
thissv <- read.event(event=caldir,t=i,var="vbsc")$vbsc
sv <- sv+thissv
}
### (8) Store the calibration factors of the current tilt: ###
cali[,k] <- truesv/sv[j,]*runs
}
### (9) Write and return the calibration data: ###
# Apply the mean of all beams if usemean==TRUE (omnidirectional sonars such as Simrad SX90):
if(usemean){
#out=c(list(cali=apply(cali,2,function(xx) rep(mean(xx),length(xx))), cal0=cali, grde=dire, ctim=unclass(Sys.time())), data[setdiff(names(data),c("dire","ftim","mtim","utim"))])
out <- c(list(cali=apply(cali,2,function(xx) rep(mean(xx),length(xx))), cal0=cali, grde=dire, ctim=unclass(Sys.time())))
}
else{
#out=c(list(cali=cali, grde=dire, ctim=unclass(Sys.time())), data[setdiff(names(data),c("dire","ftim","mtim","utim"))])
out <- c(list(cali=cali, grde=dire, ctim=unclass(Sys.time())))
}
# Write the calibration data:
if(length(directory)>0 && !identical(directory,FALSE)){
con <- file.path(directory, paste("CalibrationTable_", toupper(esnm), "_Cruise_", cruise, if(nchar(nameadd)>0) "_", nameadd, "_sv.beams", sep=""))
write.TSD(out, con=con, numt=1)
}
# Return the calibration data:
c(out, list(truesv=truesv, caldir=caldir))
}
else{
### (5) Generate the target positions and write the school file: ###
region <- get.specs.esnm(data, esnm=esnm, add=add)$region
#dr=data$asps*data$sint/2
r0 <- (j-2)*data$rres
r1 <- j*data$rres
# Total volume in which the target positions are drawn. Also see appendix 2 of the documentation of echoIBM:
V <- vol.sph(r=c(r0,r1),theta=region[,1],phi=region[,2])$volx
# The true volume backscatter, equal to N*sigma_bs/V = N*A_0/V for spherical targets:
truesv <- sgbs*N/V
# Generate the target positions and write to the school file at each run:
cat("Write the dynamic school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_dynamic.school",sep="")
for(i in seq_len(runs)){
cat("Run ",i,":\t",sep="")
### Generate positions corresponding to uniformly distributed objects in space, from which to simulate observations from the MS70 sonar. See appendix 2 of the documentation of echoIBM: ###
xyz <- runif.sph(N,r=c(r0,r1),theta=region[,1],phi=region[,2],sph.out=FALSE)
#plot3d(xyz,size=0.3)
# Add transducer position and cut the positions above the sea level:
xyz[,3] <- xyz[,3]+data$psze[1]
xyz <- xyz[xyz[,3]<0,]
# Update 'N':
thisN <- nrow(xyz)
cat(thisN,"point sources drawn\n")
# Add the target information (all rotations left out because of points targets):
schooldynamic <- list(psxf=xyz[,1],psyf=xyz[,2],pszf=xyz[,3],utim=i-1)
# Write the dynamic school data to file:
if(i==1){
write.TSD(schooldynamic, schooldynamicFile, reserve=runs-1, numt=1)
}
else{
write.TSD(schooldynamic, schooldynamicFile, numt=1, append=TRUE)
}
}
### (6) Run the simulation: ###
echoIBM(caldir, t="all", TVG.exp=2, calibrate=FALSE, noise="", max.memory=max.memory, max.radius=max.radius, method=method, esnm=esnm, parlist=list(seed=seed), keep.temp=FALSE)
### (7) Read the simulated data: ###
sv <- zeros(max(data$lenb),length(data$lenb))
for(i in seq_len(runs)){
thissv <- read.event(event=caldir,t=i,var="vbsc")$vbsc
sv <- sv+thissv
}
### (8) Write the calibration data if required: ###
calfact <- truesv/sv[j,]*runs
### (9) Write and return the calibration data: ###
# Apply the mean of all beams if usemean==TRUE (omnidirectional sonars such as Simrad SX90):
if(usemean){
#out=c(list(cali=rep(mean(calfact),length(calfact)), cal0=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
out <- c(list(cali=rep(mean(calfact),length(calfact)), cal0=calfact, ctim=unclass(Sys.time())))
}
else{
#out=c(list(cali=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
out <- c(list(cali=calfact, ctim=unclass(Sys.time())))
}
# Write the calibration data:
if(length(directory)>0 && !identical(directory,FALSE)){
con <- file.path(directory,paste("CalibrationTable_",toupper(esnm),"_Cruise_",cruise,if(nchar(nameadd)>0) "_",nameadd,"_sv.beams",sep=""))
write.TSD(out, con=con, numt=1)
}
# Return the calibration data:
c(out, list(truesv=truesv, caldir=caldir))
}
}
# The calibration procedure for target strength 'TS' is as follows: Place one target on the acoustic axis of each beam, and simulate the TS using type="ts" in echoIBM(). This includes 40*log*r for the TVG function instead of 20*log*r as for the sv. Compare the simulated 'sigma_bs' to the true 'sigma_bs' by defining calibration factors equal to true 'sigma_bs' / colSums(simulated 'sigma_bs').
else if(strff("t",type[1])){
### (1) Create the directory of the calibration simulation: ###
if(length(directory)>0 && !identical(directory,FALSE)){
calname <- paste("Calibration",esnm,"Cruise",cruise,"sigma_bs",sep="_")
}
else{
calname <- paste("Calibration",esnm,"Cruise",cruise,"sigma_bs_TEST",sep="_")
}
caldir <- paste("~/Data/echoIBM/Calibration/Events/",calname,"/",esnm,"/tsd",sep="")
if(is.na(file.info(caldir)$isdir)){
dir.create(caldir,recursive=TRUE)
}
else{
ans <- readline("The calibration case alreaddy exists. Replace? (y/n)")
if(!identical(ans,"y")){
stop("Calibration terminated")
}
else{
tempdir <- paste("~/Data/echoIBM/Calibration/Events/",calname,"/",esnm,"/temp",sep="")
if(isTRUE(file.info(tempdir)$isdir)){
unlink(tempdir, recursive = TRUE)
}
warning("The existing calibration case replaced")
}
}
### (2) Write the beam-configuration data: ###
cat("Write the beam-configuration data...\n")
write.TSD(data, paste(caldir,"/",calname,".beams",sep=""), numt=1)
### (3) Write the vessel data: ###
cat("Write the vessel data...\n")
vessel <- list(psxv=zeros(Ni),psyv=zeros(Ni),pszv=zeros(Ni),rtzv=zeros(Ni),utim=sequence(Ni)-1)
write.TSD(vessel,paste(caldir,"/",calname,".vessel",sep=""),numt=Ni)
### (4) Write the static school data: ###
schoolstatic <- list(pbpf="ps",sgbs=sgbs,obln=1,gamw=0,gaml=0)
cat("Write the static school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_static.school",sep="")
write.TSD(schoolstatic, schooldynamicFile, numt=1)
# If 'dire' is given, run through its elements and calibrate for each value, equal for all beams:
if(length(dire)>0){
cali <- NAs(Ni, length(dire))
for(k in seq_along(dire)){
cat("Calibration for tilt angle ",k," (of",length(dire),") ...\n")
# Use the input elevation angle:
thisdire <- list(dire=rep(dire[k],length(data$lenb)))
### (5) Generate the target positions and write the school file: ###
rthph <- cbind((j-1)*data$rres, data$dira, thisdire$dire)
xyz <- sph2car(rthph)
xyz[,3] <- xyz[,3]+data$psze[1]
schooldynamic <- list(psxf=xyz[,1],psyf=xyz[,2],pszf=xyz[,3],utim=sequence(Ni)-1)
cat("Write the dynamic school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_dynamic.school",sep="")
write.TSD(schooldynamic, schooldynamicFile, numt=Ni)
### (6) Run the simulation: ###
echoIBM(caldir, adds=thisdire, t="all", TVG.exp=4, calibrate=FALSE, noise="", max.memory=max.memory, method=method, max.radius=max.radius, esnm=esnm, parlist=list(seed=seed), keep.temp=FALSE)
### (7) Read the simulated data: ###
thiscali <- zeros(Ni)
for(i in seq_len(Ni)){
thissv <- read.event(event=caldir,t=i,var="vbsc")$vbsc
thiscali[i] <- sum(thissv[,i])
}
### (8) Store the calibration factors of the current tilt: ###
cali[,k] <- sgbs/thiscali
}
### (9) Write and return the calibration data: ###
# Apply the mean of all beams if usemean==TRUE (omnidirectional sonars such as Simrad SX90):
if(usemean){
#out=c(list(cali=apply(cali,2,function(xx) rep(mean(xx),length(xx))), cal0=cali, grde=dire, ctim=unclass(Sys.time())), data[setdiff(names(data),c("dire","ftim","mtim","utim"))])
out <- c(list(cali=apply(cali,2,function(xx) rep(mean(xx),length(xx))), cal0=cali, grde=dire, ctim=unclass(Sys.time())))
}
else{
#out=c(list(cali=cali, grde=dire, ctim=unclass(Sys.time())), data[setdiff(names(data),c("dire","ftim","mtim","utim"))])
out <- c(list(cali=cali, grde=dire, ctim=unclass(Sys.time())))
}
# Write the calibration data:
if(length(directory)>0 && !identical(directory,FALSE)){
con <- file.path(directory,paste("CalibrationTable_",toupper(esnm),"_Cruise_",cruise,if(nchar(nameadd)>0) "_",nameadd,"_sigma_bs.beams",sep=""))
write.TSD(out, con=con, numt=1)
}
# Return the calibration data:
c(out, list(truesv=truesv, caldir=caldir))
}
else{
### (5) Generate the target positions and write the school file: ###
rthph <- cbind((j-1)*data$rres, data$dira, data$dire)
xyz <- sph2car(rthph)
xyz[,3] <- xyz[,3]+data$psze[1]
schooldynamic <- list(psxf=xyz[,1],psyf=xyz[,2],pszf=xyz[,3],utim=sequence(Ni)-1)
cat("Write the dynamic school data...\n")
schooldynamicFile <- paste(caldir,"/",calname,"_pointTargets_dynamic.school",sep="")
write.TSD(schooldynamic, schooldynamicFile, numt=Ni)
### (6) Run the simulation: ###
echoIBM(caldir, t="all", TVG.exp=4, calibrate=FALSE, noise="", max.memory=max.memory, method=method, max.radius=max.radius, esnm=esnm, parlist=list(seed=seed), keep.temp=FALSE, cores=cores)
### (7) Read the simulated data: ###
cali <- zeros(Ni)
for(i in seq_len(Ni)){
thissv <- read.event(event=caldir,t=i,var="vbsc")$vbsc
cali[i] <- sum(thissv[,i])
}
### (8) Write the calibration data if required: ###
calfact <- sgbs/cali
# Apply the mean of all beams if usemean==TRUE (omnidirectional sonars such as Simrad SX90):
if(usemean){
#out=c(list(cali=rep(mean(calfact),length(calfact)), cal0=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
out <- c(list(cali=rep(mean(calfact),length(calfact)), cal0=calfact, ctim=unclass(Sys.time())))
}
else{
#out=c(list(cali=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
#out=c(list(cali=calfact, ctim=unclass(Sys.time())), data[setdiff(names(data),c("ftim","mtim","utim"))])
out <- c(list(cali=calfact, ctim=unclass(Sys.time())))
}
# Write the calibration data:
if(length(directory)>0 && !identical(directory,FALSE)){
con <- file.path(directory,paste("CalibrationTable_",toupper(esnm),"_Cruise_",cruise,if(nchar(nameadd)>0) "_",nameadd,"_sigma_bs.beams",sep=""))
write.TSD(out, con=con, numt=1)
}
# Return the calibration data:
c(out, list(truesv=truesv, caldir=caldir))
}
}
# Else return an error if the calibration type is not recognized:
else{
stop("Invalid calibration type (must be one of \"sv\" or \"ts\")")
}
##################################################
##################################################
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.