R/echoIBM.setVessel.R

In arnejohannesholmin/echoIBM: Simulation of multibeam sonar data

#*********************************************
#*********************************************
#' Generates and writes vessel files for all acoustic instruments of a simulation event. The function \code{echoIBM.getVesselPath} generates the vessel information for a specific esnm.
#'
#' The vessel information can be given in 4 ways: 
#' \itemize{
#'	\item nodesEarth + speed
#'	\item origin + nodesLocal + speed
#'	\item origin + heading + distance + speed
#'	\item origin + heading + duration + speed
#'	}
#'
#' @param event			A list of the following elements: (1) 'path', giving the paths to the sub-events, (2) 'esnm', giving the names of the acoustic instruments in the events (same length as 'path'), and (3) 'name', giving the name of the event.
#' @param starttime		The start time of the event, i.e., the time of the first ping/vessel position, either given as a time object or a string such as "2015-01-01 00:00:00". 
#' @param utim			The UNIX time points of the event. If not given, the starttime and pingduration will define the UNIX time points.
#' @param origin		The origin of the event, i.e., a vector of two elements giving the longitude and latitude of the first ping.
#' @param heading		The heading of the vessel in radians counter clockwise from East, either given as a single numeric, or as a vector of headings associated with vessel track segments.
#' @param distance		The distance associated with each vessel track segment.
#' @param duration		The duration of each vessel track segment in seconds.
#' @param speed			The speed of the vessel in knots, associated with each vessel track segment.
#' @param nodesLocal	The nodes defining the vessel track segments, given as local nodes relative to the origin (can include the origin c(0, 0))
#' @param nodesEarth	The nodes defining the vessel track segments, given as global nodes in the coordinate system of the earth relative to the origin (can include the origin c(0, 0)).
#' @param pingduration	The duration of the pings in seconds. Currently only one single fixed value is allowed.
#' @param heave			The heave of the vessel in meters, either given as a single value, a vector of length equal to the number of vessel track segments, or alternatively a funciton of the number of time steps such as rnorm() or the more appropriate function(x) {set.seed(x); runif(x)}, which sets the seed as the number of time steps.
#' @param rtxv			The pitch of the vessel in radians positive for uppwards pitch (bow lifting). See info.TSD("rtxv")
#' @param rtyv			The roll of the vessel in radians positive for starboard side tilted down. See info.TSD("rtyv")
#' @param rtzv			The orientation of the vessel in radians counter clockwise, IN ADDITION to the heading of the vessel specified in \code{heading}. See info.TSD("rtzv")
#' @param ...			Data overriding the variables generated by the funciton.
#'
#' @return
#'
#' @examples
#' \dontrun{}
#'
#' @importFrom TSD write.TSD
#'
#' @export
#' @rdname echoIBM.setup
#'
echoIBM.setVessel <- function(
	event, 
	starttime = "2015-01-01 00:00:00", 
	utim = NULL, 
	origin = c(0, 0), 
	heading = 0, # East
	distance = NULL, 
	duration = 1, 
	speed = NULL, 
	nodesLocal = NULL, 
	nodesEarth = NULL, 
	pingduration = 1,
	numt = NULL, 
	heave = 0, 
	rtxv = 0, 
	rtyv = 0, 
	rtzv = 0, 
	offset = NULL, 
	freeze = NULL, 
	...){
	
	############### LOG: ###############
	# Start: 2017-03-29 - Clean version.
	
	# Save the input variables:
	dotList <- list(...)
	
	vessel <- echoIBM.getVesselPath(starttime=starttime, origin=origin, heading=heading, distance=distance, speed=speed, duration=duration, nodesLocal=nodesLocal, nodesEarth=nodesEarth, utim=utim, pingduration=pingduration, numt=numt, heave=heave, rtxv=rtxv, rtyv=rtyv, rtzv=rtzv)
	
	# Add data:
	vessel <- replaceKeepDim(vessel, dotList, esnm="")
	
	# Expand all variables to equal length for faster reading and writing:
	vessel <- lapply(vessel, rep, length.out=max(sapply(vessel, length)))
	
	# Define files: 
	vesselfiles <- file.path(event$path, paste0(event$name, ".vessel"))
	names(vesselfiles) <- event$esnm
	
	# Repeat the vessel data:
	vesselAll <- rep(list(vessel), length(vesselfiles))
	names(vesselAll) <- event$esnm
	# Apply offset on the vessel data:
	nameMatch <- intersect(names(offset), names(vesselAll))
	for(name in nameMatch){
		vesselAll[[name]] <- addVesselOffset(vesselAll[[name]], offset=offset[[name]])
	}
	# Freeze the vessel data:
	nameMatch <- intersect(names(freeze), names(vesselAll))
	for(name in nameMatch){
		vesselAll[[name]] <- freezeVesselPosition(vesselAll[[name]], freeze=freeze[[name]])
	}
	
	# Run through the events and write the vessel file:
	lapply(vesselfiles, function(x) suppressWarnings(dir.create(dirname(x), recursive=TRUE)))
	#lapply(vesselfiles, function(file) write.TSD(vessel, con=file, numt=vessel$numt[1]))
	lapply(seq_along(vesselfiles), function(i) write.TSD(vesselAll[[i]], con=vesselfiles[i], numt=vesselAll[[i]]$numt[1]))
	
	# Add the files to the output:
	return(vesselfiles)
}
#'
#' @importFrom TSD global2car car2global ftim2utim
#'
#' @export
#' @rdname echoIBM.setup
#'
echoIBM.getVesselPath <- function(
	starttime = "2015-01-01 00:00:00", 
	origin = NULL, 
	heading = NULL, 
	distance = NULL, 
	speed = NULL, 
	duration = NULL, 
	nodesLocal = NULL, 
	nodesEarth = NULL, 
	utim = NULL, 
	pingduration = 1,
	numt = NULL, 
	heave = 0, 
	rtxv = 0, 
	rtyv = 0, 
	rtzv = 0){
		
	############### LOG: ###############
	# Start: 2017-03-29 - Clean version.
	# Function for applying randmoness to vessel dynamics:
	setVesselDynamics <- function(x, x0, nums, segmentind){
		# Set the default values, usually 0 but for rtzv it should be the heading:
		x0 <- rep(x0, length.out=nums)[segmentind]
		# Apply a function of the number of time steps, or simply repeat the numeric value, which is assumed to correspond to the vessel track segments:
		if(is.numeric(x)){
			x <- rep(x, length.out=nums)[segmentind]
		}
		else if(is.function(x)){
			x <- x(length(segmentind))
		}
		# Add the default and the altered values:
		x + x0
	}
	
	# Set the speed if missing:
	if(length(speed)==0){
		if(length(numt)==0){
			warning("Speed not given and was set to 10 knots = 10 * 1852 / 3600 = 5.14 m/s")
		}
	}
	# Convert from knots to meters per second:
	speed <- speed * 1852 / 3600
	
	# If given as nodes in the coordinate system of the Earth, convert to local nodes:
	if(length(nodesEarth)){
		if(length(origin)==0){
			origin <- nodesEarth[1, 1:2]
		}
		nodesLocal <- global2car(nodesEarth, origin=origin)
	}
	
	# If local nodes and origin are given or deduced, extract heading and distance
	if(length(nodesLocal)){
		if(length(origin)==0){
			origin <- c(0, 0)
		}
		segments <- diff(nodesLocal[,1:2])
		heading <- atan2(segments[,2], segments[,1])
		distance <- sqrt(rowSums(segments^2))
		duration <- distance / speed
	}
	
	### # Add the origin to the nodesLocal if not present:
	### if(length(nodesLocal) && !all(nodesLocal[1, 1:2] == c(0,0))){
	### 	nodesLocal <- rbind(c(0,0), nodesLocal)
	### }
	### # If local nodes and origin are given or deduced, extract heading and distance
	### if(length(nodesLocal) && length(origin)){
	### 	segments <- diff(nodesLocal[,1:2])
	### 	heading <- atan2(segments[,2], segments[,1])
	### 	distance <- sqrt(rowSums(segments^2))
	### 	duration <- distance / speed
	### }
	# If heading, duration and origin is not given, issue an error:
	if(!any(length(duration), length(distance)) && !all(length(origin), length(heading))){
		stop("Dynamics need to be given either as (1) 'nodesEarth', (2) 'nodesLocal' and 'origin', or (3) 'origin', 'heading', and one of 'duration' and 'distance'")
	}
	# If these are present, extract the local nodes:
	else if(length(nodesLocal)==0){
		if(length(distance)==0){
			distance <- duration * speed
		}
		nodesLocal <- cbind(c(0, cumsum(distance * cos(heading))), c(0, cumsum(distance * sin(heading))))
	}
	
	# If origin, heading, speed and distance is given, extract duration instead of distance:
	if(length(numt)){
		totalDistance <- sum(distance)
		totalDuration <- numt * pingduration
		speed <- totalDistance / totalDuration
		duration <- distance / speed
	}
	else if(length(duration)==0 && length(origin) && length(heading) && length(speed) && length(distance)){
		duration <- distance / speed
	}
	
	totalDuration <- sum(duration)
	
	
	# Get numt, and the times as spread evenly along the time span:
	if(length(utim)==0){
		utim0 <- ftim2utim(starttime)
		if(length(numt)){
			utim <- seq(0, totalDuration, length.out=numt + 1)
		}
		else{
			utim <- seq(0, totalDuration, pingduration)
			numt <- length(utim) - 1
		}
		utim <- utim[-1] - diff(utim)/2
	
	}
	else{
		numt <- length(utim)
		utim0 <- utim[1]
		utim <- utim - utim0
	}
	
	
	nums <- length(duration)
	
	
	# Get cummulative duration:
	cduration <- cumsum(duration)
	# And intervals of duration:
	cduration0 <- c(0, cduration)
	
	# Get the indices of the segments of each time step:
	segmentind <- findInterval(utim, cduration0, rightmost.closed=TRUE)
	# And the duration in the present segment for each time step:
	durationInSegment <- utim - cduration0[segmentind]
	# Add the first time to the output time:
	utim <- utim0 + utim
	
	# repeat speed, heave, rtxv and rtyv by segments:
	speed <- rep(speed, nums)[segmentind]
	ispv <- speed[segmentind]

	# Set vessel dynamics not derived elsewhere in the function:
	heave <- setVesselDynamics(heave, 0, nums, segmentind)
	rtxv <- setVesselDynamics(rtxv, 0, nums, segmentind)
	rtyv <- setVesselDynamics(rtyv, 0, nums, segmentind)
	rtzv <- setVesselDynamics(rtzv, heading, nums, segmentind)

	
	# Get the positions
	psxy <- nodesLocal[segmentind,1:2] + durationInSegment * speed[segmentind] * cbind(cos(heading[segmentind]), sin(heading[segmentind]))
	lonlat <- car2global(psxy, origin=origin)
	
	out <- list(ispv=ispv, psxv=psxy[,1], psyv=psxy[,2], pszv=heave, rtxv=rtxv, rtyv=rtyv, rtzv=rtzv, lonv=lonlat[,1], latv=lonlat[,2])
	# Get sailed distance in nautical miles:
	nmi <- 1852
	sadv <- sqrt(diff(out$psxv)^2 + diff(out$psyv)^2) / nmi
	out$sadv <- c(0, cumsum(sadv))
	
	out <- c(list(utim=utim, numt=numt), out, list(lon0=origin[1], lat0=origin[2]))
	
	# Repeat to maximum dimension:
	out <- repToExtreme(out, ndim=1, skip.len=1)
	
	return(out)
}
addVesselOffset <- function(vessel, offset){
	for(name in names(offset)){
		vessel[[name]] <- vessel[[name]] + offset[[name]]
	}
	
	# Update geographical positions:
	lonlat <- car2global(cbind(vessel$psxv, vessel$psyv), origin=c(vessel$lon0[1], vessel$lat0[1]))
	vessel$lonv <- lonlat[,1]
	vessel$latv <- lonlat[,2]
	
	vessel
}
freezeVesselPosition <- function(vessel, freeze){
	for(name in c("psxv", "psyv")){
		# Repeat the specified time step:
		numt <- length(vessel[[name]])
		if(freeze>0 && freeze<1){
			freeze <- freeze * numt
		}
		freeze <- round(freeze)
		vessel[[name]] <- rep(vessel[[name]][freeze], length.out=numt)
	}
	
	# Update geographical positions:
	lonlat <- car2global(cbind(vessel$psxv, vessel$psyv), origin=c(vessel$lon0[1], vessel$lat0[1]))
	vessel$lonv <- lonlat[,1]
	vessel$latv <- lonlat[,2]
	
	vessel
}