R/bottom.contact.r
In jae0/netmensuration: Package that interacts with bio and associated projects to handle net mensuration data streams
Defines functions
bottom.contact
Documented in
bottom.contact
#' @title bottom.contact
#' @description Calculates indicies of touchdown and liftoff
#' @import tcltk geosphere lubridate
#' @return dataframe with various results
#' @export
bottom.contact = function( x, bcp, debugrun=FALSE ) {
require(tcltk)
# all timestamps must be in Posix/UTC
# range comparisons seem to fail when they are not
tzone = "UTC"
if (debugrun) {
debug.plot = TRUE
browser()
}
debug.plot = FALSE
n.min.required = 30
nx = nrow(x)
O = list() # output list
O$id = bcp$id
O$error.flag = NA
O$good = rep(TRUE, nx) # rows that will contain data that passes each step of data quality checks
if (length (which( (!is.na(x$depth)))) < n.min.required ) return( NULL )
##--------------------------------
# sort in case time is not in sequence
# timestamps have frequencies higher than 1 sec .. duplciates are created and this can pose a problem
x = x[order( x$timestamp ) ,]
x$ts = as.numeric( difftime( x$timestamp, min(x$timestamp), units="secs" ) )
O$plotdata = x # save incoming data after creation of time stamps and ordering
## PRE-FILTER 1
if ( exists( "double.depth.sensors", bcp) ) {
# two depth sensors were used simultaneously but they are not calibrated!
# remarkably hard to filter this out with any reliability at a later stage ..
om = modes( x$depth )
oaoi = range( which( x$depth <= om$ub2 & x$depth >= om$lb2))
oj = oaoi[1]:oaoi[2]
oi = interpolate.xy.robust( x[oj, c("ts", "depth")], method="loess", trim=0.05, probs=c(0.05, 0.95), target.r2=0.5 )
oo = which( (x$depth[oj] - oi) < 0)
x$depth[oj[oo]] = NA
if (any( is.na( x$depth))) x$depth = interpolate.xy.robust( x[, c("ts", "depth")], method="simple.linear" )
}
## PRE-FILTER 2
if (all(!is.na(bcp$time.gate))) {
# simple time-based gating if a time range is provided..
bcp$trange.max = (bcp$tdif.max+5)
O$good = bottom.contact.gating.time ( Zt=x$timestamp, good=O$good, bcp=bcp )
x$depth[ !O$good ] = NA
}
## PRE-FILTER 3
# simple depth-based gating
if( !any(x$depth > bcp$depth.min,na.rm=T)) return( NULL )
O$good = bottom.contact.gating.depth ( Z=x$depth, good=O$good, bcp=bcp)
x$depth[ !O$good ] = NA
## PRE-FILTER 4
# time and depth-based gating
mm = modes( x$depth ) # naive first estimate of location of most data depths
if (mm$lb == mm$ub || is.na(mm$sd)) return(NULL) # there is no depth variation in the data ... likely a bad data series
mm.i = which( x$depth > (mm$lb2+bcp$depth.range[1]) & x$depth < (mm$ub2 + bcp$depth.range[2]) )
O$good[ setdiff(1:nx, mm.i)] = FALSE
O$good[mm.i] = TRUE
## RANGE CHECKS
x.timerange = range( x$timestamp[O$good], na.rm=TRUE )
x.dt = difftime( x.timerange[2], x.timerange[1], units="mins" )
if ( x.dt < (bcp$tdif.min ) ) {
O$error.flag = "Time track is too short?"
return(O)
}
if ( x.dt > (bcp$tdif.max+10) ) { # +10 = 5 min on each side
# data vector is too long ... truncate
mm.r = range( mm.i )
mm.dt = difftime( x$timestamp[mm.r[2]], x$timestamp[mm.r[1]], units="mins" )
if ( mm.dt > (bcp$tdif.max+10) ) {
O$error.flag = "Time track is too long?"
}
}
if ( sd(x$depth, na.rm=TRUE) < bcp$eps.depth ) {
O$error.flag = "not enough variability in data?"
return(O)
}
## inSANITY CHECKS
# sometimes multiple tows exist in one track ...
# over-smooth depth to capture strange tows
x$dsm = interpolate.xy.robust( x[, c("ts", "depth")], method="sequential.linear", trim=0.05 )
x$dsm = interpolate.xy.robust( x[, c("ts", "dsm")], method="local.variance", trim=0.05 )
x$dsm = interpolate.xy.robust( x[, c("ts", "dsm")], method="sequential.linear", trim=0.05 )
x$dsm = interpolate.xy.robust( x[, c("ts", "dsm")], method="loess", trim=0.05 )
zz = rep( 0, nx )
zz[ which( x$dsm < ( mm$mode / 3 )) ] = 1 # flag shallow areas
inc.depth = abs( diff( x$dsm ) )
# also capture strong noise - very obviously wrong data
rapid.depth.changes = which( inc.depth > bcp$maxdepthchange )
if ( length( rapid.depth.changes ) > 0 ) {
zz[ rapid.depth.changes ] = 1
zz[ rapid.depth.changes-1 ] = 1 # include adjecent points to remove
zz[ rapid.depth.changes+1 ] = 1
}
dzz = diff(zz)
bnds = c(1, which( dzz != 0 ), nx )
if ( length (bnds) > 2 ) {
# i.e. , contaminated by noise or multiple tows
segs = diff(bnds) # length of each segment
longest = which.max( segs )
gg = bnds[longest]:bnds[(longest+1)]
bad = setdiff( 1:nx, gg)
O$good[bad] = FALSE
}
# fix tails
mm = modes( x$depth[ O$good ] ) # recompute after gating
mm.bot = which ( x$depth >= mm$lb2 & x$depth <= mm$ub2 ) # first est of bottom
if(length(mm.bot) == 0) return(NULL) # not enough variability in data
mm.med = floor( median( mm.bot, na.rm=TRUE ) )
# fix tails when the tails are not a single smooth tail:
llim = mm$mode + bcp$depth.range[1]
for ( ll in mm.med:1 ) {
if ( is.finite( x$depth[ll] ) ) {
if ( x$depth[ll] < llim ) break()
}
}
for ( uu in mm.med:nrow(x)) {
if ( is.finite( x$depth[uu] ) ) {
if ( x$depth[uu] < llim ) break()
}
}
todrop = c( 1:ll, uu:nrow(x) )
x$depth[ todrop ] = NA
O$good[ todrop] = FALSE
if ( sd(x$depth, na.rm=TRUE) < bcp$eps.depth ) return(NULL) # not enough variability in data
## ------------------------------
## MAIN NOISE/INTERPOLATION FILTER
# Some filtering of noise from data and further focus upon area of interest based upon time and depth if possible
res = NULL
res = try( bottom.contact.filter.noise ( x=x, good=O$good, bcp=bcp, debug=debugrun ), silent =TRUE )
if ( "try-error" %in% class(res) || is.null(res$depth.smoothed) ) {
x$depth = jitter( x$depth )
res = try( bottom.contact.filter.noise ( x=x, good=O$good, bcp=bcp, debug=debugrun ), silent =TRUE )
}
if ( !"try-error" %in% class(res) ) {
if (!is.null(res$depth.smoothed)) {
if ( cor( x$depth, res$depth.smoothed, use="pairwise.complete.obs") > 0.999 || is.na(cor( x$depth, res$depth.smoothed, use="pairwise.complete.obs") )) {
bcp$noisefilter.target.r2 = bcp$noisefilter.target.r2 - 0.1
x$depth = jitter( x$depth )
res = try( bottom.contact.filter.noise ( x=x, good=O$good, bcp=bcp, debug=debugrun ), silent =TRUE )
}}}
x$depth.smoothed = x$depth
if ( ! "try-error" %in% class( res) ) {
# retain good/bad information only for the interval where there is data
# .. this effectively does the filtering in the area of interest only and retains the tails
# for further influence later
ig = range( which( res$good ))
igg = ig[1]:ig[2]
O$good[igg] = res$good[igg]
x$depth[ !O$good ] = NA
if (!is.null(res$depth.smoothed)) x$depth.smoothed= res$depth.smoothed
}
if(sum(x$depth-min(x$depth,na.rm=T),na.rm=T)==0) return( NULL )
if(sum(O$good)==0) return( NULL )
if (sd(x$depth, na.rm=TRUE) < bcp$eps.depth ) return(NULL) # not enough variability in data
if (sd(x$depth.smoothed, na.rm=TRUE) < bcp$eps.depth ) return(NULL) # not enough variability in data
x.timerange = range( x$timestamp[O$good], na.rm=TRUE )
x.dt = difftime( x.timerange[2], x.timerange[1], units="mins" )
if ( x.dt < ( bcp$tdif.min ) ) {
O$error.flag = "Not enough data?"
return(O)
}
if ( x.dt > (bcp$tdif.max + 10 ) ) {
O$error.flag = "Too much data?"
return(O)
}
## FINAL DATA GATING
# variance gating attempt after the data has been cleaned as much as possible
O$variance.method0 = NA
O$variance.method1 = NA
O$variance.method.indices = NA
res = NULL
res = try( bottom.contact.gating.variance ( x, O$good, bcp ), silent =TRUE )
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1 ) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$variance.method0 = res$bc0
O$variance.method1 = res$bc1
O$variance.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 )
bad = which( x$timestamp < res$bc0 | x$timestamp > res$bc1 )
if (length( bad) > 0) O$good[ bad ] = FALSE
x$depth[ !O$good ] = NA
} }
}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "gray"
points( depth~ts, x[ O$variance.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## NOTE::: From this point on, O$good is now complete
## -- it contains indices of time and depth-based gating as well as the variance based gating
# finalize selection of area of interest (based upon gating, above)
aoi.range = range( which( O$good ) )
aoi.mid = floor( mean( aoi.range ) ) # approximate midpoint
aoi.min = aoi.range[1]
aoi.max = aoi.range[2]
O$aoi = aoi.min:aoi.max # stored for use in the linear method where we need to recover the left-most index
##--------------------------------
# Modal method: gating by looking for modal distribution and estimating sd of the modal group in the data
# first by removing small densities ( 1/(length(i)/nb) ) and by varying the number of breaks in the histogram
# until a target number of breaks, nbins with valid data are found
O$modal.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$modal.method1 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$modal.method.indices = NA
sm0 = x[ O$aoi, c("depth.smoothed", "timestamp", "ts" ) ] # send filtered data ... continuity not important .. order is important
res = NULL
res = try( bottom.contact.modal( sm=sm0, bcp ), silent=TRUE )
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1 ) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$modal.method0 = res$bc0 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$modal.method1 = res$bc1 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$modal.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 ) # x correct
} }
}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "green"
points( depth~ts, x[ O$modal.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## ----------------------------
## Smooth method: using smoothed data (slopes are too unstable with raw data),
## compute first derivatives to determine when the slopes inflect
O$smooth.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$smooth.method1 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$smooth.method.indices = NA
sm0 = x[ O$aoi, c("depth.smoothed", "timestamp", "ts")] # Send all data within the aoi --- check this .. order is important
res = NULL
res = try(
bottom.contact.smooth( sm=sm0, bcp=bcp ) , silent =TRUE)
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$smooth.method0 = res$bc0 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$smooth.method1 = res$bc1 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$smooth.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 ) # x correct
} }
}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "blue"
points( depth~ts, x[ O$smooth.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## ----------------------------
## maxdepth method: looking for the max depth near the areas of interest (left, right)
O$maxdepth.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$maxdepth.method1 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$maxdepth.method.indices = NA
sm0 = x[, c("depth.smoothed", "timestamp", "ts")] # Send all data within the aoi --- check this .. order is important
sm0$depth[ !O$good ] = NA
sm0 = sm0[ O$aoi, ]
bcmethods=c( "smooth.method", "modal.method" )
res = NULL
res = try( bottom.contact.maxdepth( sm=sm0, O=O, bcmethods=bcmethods, bcp=bcp ) , silent=TRUE )
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if(length(res$bc1)>0 && length(res$bc0)>0){
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1 ) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$maxdepth.method0 = res$bc0 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$maxdepth.method1 = res$bc1 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$maxdepth.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 ) # x correct
}}
}}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "yellow"
points( depth~ts, x[ O$maxdepth.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## ---------------------------
## Linear method: looking at the intersection of three lines (up, bot and down)
O$linear.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$linear.method1 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$linear.method.indices = NA
sm0 = x[, c("depth.smoothed", "timestamp", "ts")] # Send all data within the aoi --- check this .. order is important
sm0$depth[ !O$good ] = NA
sm0 = sm0[ O$aoi, ]
bcmethods=c( "smooth.method", "modal.method", "linear.method" )
res = NULL
res = try( bottom.contact.linear( sm=sm0, O=O, bcmethods=bcmethods, bcp=bcp ) , silent=TRUE )
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$linear.method0 = res$bc0 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$linear.method1 = res$bc1 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$linear.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 ) # x correct
} }
}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "red"
points( depth~ts, x[ O$linear.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## ---------------------------
O$manual.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$manual.method1 = NA ### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
if ( bcp$user.interaction ) {
satisified = F
while(!satisified){
#open plot window, user can modify the numbers to fit their screen
x11(width = 25, height = 15)
#Define Margins for multiple y-axis
par(mar=c(5, 8, 4, 4) + 0.1)
grange = min(which(O$good)):max(which(O$good))
if("wingspread" %in% names(O$plotdata)){
#Plot the door spread
ylim=c(min(O$plotdata$wingspread[grange], na.rm = TRUE)-1,max(O$plotdata$wingspread[grange], na.rm = TRUE)+1)
if(is.finite(ylim)){
plot(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$wingspread[grange]))], O$plotdata$wingspread[grange][which(!is.na(O$plotdata$wingspread[grange]))], axes=F, ylim=ylim, xlab="", ylab="",type="p",col="#0000FF1A", main="",xlim=c(min(O$plotdata$timestamp[grange]), max(O$plotdata$timestamp[grange])))
if(length(unique(na.omit(O$plotdata$wingspread[grange]))) > 5 ){
smoothingSpline = smooth.spline(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$wingspread[grange]))], O$plotdata$wingspread[grange][which(!is.na(O$plotdata$wingspread[grange]))], spar=.5)
lines(smoothingSpline, col="blue")
}
abline( h = c(median(O$plotdata$wingspread[grange], na.rm = TRUE)), col = "blue", lty =2 )
axis(2,col="blue", col.lab = "blue", col.axis = "blue", lwd=1)
mtext(2,text="Wing Spread",col = "blue",line=2)
}
}
yl=c(min(O$plotdata$opening[grange], na.rm = TRUE)-1,max(O$plotdata$opening[grange], na.rm = TRUE)+1)
if(!is.numeric(yl[1]) | is.na(yl[1]) ) yl[1] = 0
if(!is.numeric(yl[2]) | is.na(yl[2]) ) yl[2] = 10
#Plot the opening
if("opening" %in% names(O$plotdata)){
if(length(which(!is.na(O$plotdata$opening[grange]))) > 10){
par(new=T)
plot(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$opening[grange]))], O$plotdata$opening[grange][which(!is.na(O$plotdata$opening[grange]))], axes=F, ylim=yl, xlab="", ylab="",type="p",col="#A52A2A1A", main="",xlim=c(min(O$plotdata$timestamp[grange]), max(O$plotdata$timestamp[grange])))
smoothingSpline = smooth.spline(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$opening[grange]))], O$plotdata$opening[grange][which(!is.na(O$plotdata$opening[grange]))], spar=.5)
lines(smoothingSpline, col="brown")
axis(2, ylim=yl,col = "brown",col.lab = "brown", col.axis = "brown",lwd=1,line=3.5)
mtext(2,text="Opening", col = "brown", line=5.5)
}
}
#Plot the depth (Clickable)
par(new=T)
plot(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$depth[grange]))], O$plotdata$depth[grange][which(!is.na(O$plotdata$depth[grange]))], axes=F, ylim=rev(c(min(O$plotdata$depth[grange], na.rm = TRUE)-1,max(O$plotdata$depth[grange], na.rm = TRUE)+1)), xlab="", ylab="", type="l", main=O$id,xlim=c(min(O$plotdata$timestamp[grange]), max(O$plotdata$timestamp[grange])),lwd=1)
axis(4, ylim=rev(c(min(O$plotdata$depth[grange], na.rm = TRUE)-1,max(O$plotdata$depth[grange], na.rm = TRUE)+1)),lwd=1)
mtext(4,text="Depth", line = 2)
#Set up time axis plotting
abli = seq(O$plotdata$timestamp[grange][1],O$plotdata$timestamp[grange][length(O$plotdata$timestamp[grange])], "1 min")
#Draw the time axis
axis.POSIXct(1, at = abli, format = "%H:%M:%S", labels = TRUE)
ablit = unique(as.character(lubridate::date(abli)))
mtext(text = ablit, side=1, col="black", line=2)
abline( v = abli, col = "lightgrey")
#Visualize neutral haulback
#Find positions of going backward toward origin. Nuetral haulback
disfromorigin = rep(0, 1, length(O$plotdata$latitude[grange]))
for(k in 1:(length(O$plotdata$latitude[grange])-1)){
p1 = c(O$plotdata$longitude[grange][1],O$plotdata$latitude[grange][1])
p2 = c(O$plotdata$longitude[grange][k+1],O$plotdata$latitude[grange][k+1])
if(!is.null(p1) & !is.null(p2)){
if(!is.na(p1) & !is.na(p2))
{
disfromorigin[k+1] = geosphere::distHaversine(p1 , p2, r=6378137)
}
}
}
neu = rep(0, 1, length(O$plotdata$latitude[grange]))
for(k in 1:(length(disfromorigin)-1)){
neu[k] = disfromorigin[k+1]-disfromorigin[k]
if(k > 1){
if(neu[k-1]<0 & neu[k]>=0) neu[k-1] = 0
}
}
neuhaul = F
indn = which(neu < 0)
if(length(indn) > 10){
print(paste("Neutral Haulback! ", O$id))
neuhaul = T
}
if(length(indn)>0){
abline(v = O$plotdata$timestamp[grange][indn], col = "red")
}
id = identify(O$plotdata$timestamp[grange], O$plotdata$depth[grange], labels = c(as.character(O$plotdata$timestamp[grange])),n = 1, pos = TRUE)
start = O$plotdata$timestamp[grange][id$ind]
points(O$plotdata$timestamp[grange][id$ind], O$plotdata$depth[grange][id$ind], col = "darkgreen", pch = 19)
#Code that allow identifying clicks
id2 = identify(O$plotdata$timestamp[grange], O$plotdata$depth[grange], labels = c(as.character(O$plotdata$timestamp[grange])),n = 1, pos = TRUE)
end = O$plotdata$timestamp[grange][id2$ind]
points(O$plotdata$timestamp[grange][id2$ind], O$plotdata$depth[grange][id2$ind], col = "darkgreen", pch = 19)
##Code taken from bottom.contact.plot.r
legendtext = NULL
legendcol = NULL
legendpch = NULL
if (all(is.finite(c( O$variance.method0, O$variance.method1) ) ) ) {
mcol = "pink"
# points( depth~ts, x[ O$variance.method.indices, ], pch=20, col=mcol, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$variance.method.indices)], col=mcol, lty="solid", lwd=1.2)
abline (v=O$plotdata$timestamp[max(O$variance.method.indices)], col=mcol, lty="solid", lwd=1.2)
DT = as.numeric( difftime( O$variance.method1, O$variance.method0, units="mins" ) )
legendtext = c( legendtext, paste( "variance gate: ", round( DT, 2), "" ) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20 )
}
if (all(is.finite( c(O$modal.method0, O$modal.method1 ) ) ) ) {
mcol = "red" # colour for plotting
# points( depth~ts, x[O$modal.method.indices,], col=mcol, pch=20, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$modal.method.indices)], col=mcol, lty="dashed")
abline (v=O$plotdata$timestamp[max(O$modal.method.indices)], col=mcol, lty="dashed")
DT = as.numeric( difftime( O$modal.method1, O$modal.method0, units="mins" ) )
legendtext = c( legendtext, paste( "modal: ", round( DT, 2) ) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20)
}
if ( all(is.finite( c( O$smooth.method0, O$smooth.method1) ) ) ) {
mcol = "blue"
# points( depth~ts, x[O$smooth.method.indices,], col=mcol, pch=20, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$smooth.method.indices)], col=mcol, lty="dotdash", lwd=1.5)
abline (v=O$plotdata$timestamp[max(O$smooth.method.indices)], col=mcol, lty="dotdash", lwd=1.5)
DT = as.numeric( difftime( O$smooth.method1, O$smooth.method0, units="mins" ) )
legendtext = c(legendtext, paste( "smooth: ", round(DT, 2)) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20)
}
if (all(is.finite( c(O$linear.method0, O$linear.method1) )) ) {
mcol ="green"
# points( depth~ts, x[O$linear.method.indices,], col=mcol, pch=20, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$linear.method.indices)], col=mcol, lty="twodash")
abline (v=O$plotdata$timestamp[max(O$linear.method.indices)], col=mcol, lty="twodash")
DT = as.numeric( difftime( O$linear.method1, O$linear.method0, units="mins" ) )
legendtext = c( legendtext, paste( "linear: ", round( DT, 2) ) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20)
}
if (all(is.finite( c( O$maxdepth.method0, O$maxdepth.method1) )) ) {
mcol ="orange"
# points( depth~ts, x[O$linear.method.indices,], col=mcol, pch=20, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$maxdepth.method.indices)], col=mcol, lty="solid")
abline (v=O$plotdata$timestamp[max(O$maxdepth.method.indices)], col=mcol, lty="solid")
DT = as.numeric( difftime( O$maxdepth.method1, O$maxdepth.method0, units="mins" ) )
legendtext = c( legendtext, paste( "maxdepth: ", round( DT, 2) ) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20)
}
if ( !( is.null( legendtext))) legend( "top", legend=legendtext, col=legendcol, pch=legendpch, bg = "white" )
#END code taken from bottom.contact
#User interface
goodans = FALSE
res <- tkmessageBox(title = "ManualTD",
message = "Are you happy with these values?", icon = "info", type = "yesno")
if(grepl("yes", res)){
satisified = TRUE
goodans = TRUE
O$manual.method0 = start
O$manual.method1 = end
O$manual.method.indices = which( O$plotdata$timestamp >= O$manual.method0 & O$plotdata$timestamp <= O$manual.method1 )
# result = rbind(result, row)
# if(!file.exists(file)){
# write.table(row, file = file, sep = ",", row.names = F, quote = F)
# }else{
# write.table(row, file = file, sep = ",", append = T, col.names = F, row.names = F, quote = F )
# }
}
if(grepl("no", res)){
satisified = FALSE
goodans = TRUE
}
if(!goodans)print("Did not receive valid command, please try again.")
dev.off()
}#END plotting loop while not satisfied
if(is.na(O$manual.method0) || is.na(O$manual.method1)){
bcp$user.interaction = FALSE
}
else{
if(is.null(bcp$station)) station = unlist(strsplit(bcp$id, "\\."))[4]
else station = bcp$station
bcp$id = bcp$trip
mf = NULL
if(!is.null(bcp$from.manual.archive)) mf = file.path(bcp$from.manual.archive, "clicktouchdown_all.csv")
if(!is.null(bcp$from.manual.file)) mf = bcp$from.manual.file
if(!is.null(mf)){
manualclick = NULL
if(file.exists(mf)){
manualclick = read.csv(mf, as.is=TRUE)
if(bcp$datasource == "lobster"){
sta.ind = which(manualclick$station == station & manualclick$trip == bcp$trip)
}
else{
sta.ind = which(manualclick$station == station & manualclick$year == bcp$YR)
}
if(length(sta.ind) == 0){
print(bcp$trip)
manualclick = rbind(manualclick, data.frame(station = station, start = unlist(as.character(O$manual.method0)), end = unlist(as.character(O$manual.method1)), depth = mean( x$depth, na.rm=TRUE ), year = bcp$YR, trip = bcp$id))
}
else{
manualclick$station[sta.ind] = station
manualclick$start[sta.ind] = as.character(O$manual.method0)
manualclick$end[sta.ind] = as.character(O$manual.method1)
manualclick$depth[sta.ind] = mean( x$depth, na.rm=TRUE )
manualclick$year[sta.ind] = bcp$YR
manualclick$trip[sta.ind] = bcp$id
}
}
else{
manualclick = data.frame(station = station, start = as.character(O$manual.method0), end = as.character(O$manual.method1), depth = mean( x$depth, na.rm=TRUE ), year = bcp$YR, trip = bcp$id)
}
write.csv(manualclick, mf, row.names = FALSE )
}
}
}
#BC- Added condition in case user interaction is desired, no need to do this step.
if (!is.null(bcp$from.manual.archive) && !bcp$user.interaction) {
print( "Loading values from previously generated .csv")
if(file.exists(file.path(bcp$from.manual.archive, "clicktouchdown_all.csv"))){
manualclick = read.csv(file.path(bcp$from.manual.archive, "clicktouchdown_all.csv"), as.is=TRUE)
station = unlist(strsplit(bcp$id, "\\."))[4]
sta.ind = which(manualclick$station == station & manualclick$year == bcp$YR)
if(length(sta.ind == 1)){
O$manual.method0 = lubridate::ymd_hms(manualclick$start[sta.ind], tz = "UTC")
O$manual.method1 = lubridate::ymd_hms(manualclick$end[sta.ind], tz = "UTC")
}
}
}
O$means0 = NA ### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$means1 = NA ### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
bcmethods = c("manual.method", "variance.method", "smooth.method", "modal.method", "maxdepth.method", "linear.method", "means" )
standard = which( bcmethods=="manual.method") # gold standard
direct = which( bcmethods %in% c("smooth.method", "modal.method", "linear.method", "maxdepth.method" ) )
imeans = which( bcmethods == "means" )
# must be careful as timestamp is being converted to text and will lose tzone ... need to reset to correct timezone:
bcm0 = paste(bcmethods, "0", sep="")
bcm1 = paste(bcmethods, "1", sep="")
# recovert to time zone of incoming data as time zone is lost with the transpose
tmp0 = lubridate::ymd_hms( t( as.data.frame(O[ bcm0 ]) ) ) # UTC
tmp1 = lubridate::ymd_hms( t( as.data.frame(O[ bcm1 ]) ) ) # UTC
bottom0.mean = mean(tmp0, na.rm=TRUE)
bottom1.mean = mean(tmp1, na.rm=TRUE)
bottom0.sd = max( 1, sd( tmp0, na.rm=TRUE ) ) # in seconds
bottom1.sd = max( 1, sd( tmp1, na.rm=TRUE ) )
dflag0 = rowSums(as.matrix(dist(tmp0[direct], upper=TRUE )), na.rm=TRUE) # which are most extreme
dflag1 = rowSums(as.matrix(dist(tmp1[direct], upper=TRUE )), na.rm=TRUE) # which are most extreme
tooextreme0 = which.max( dflag0 )
tooextreme1 = which.max( dflag1 )
trimmed0 = trimmed1 = direct
if (length( which(dflag0 > 0) ) > 1 ) trimmed0 = trimmed0[-tooextreme0 ]
if (length( which(dflag1 > 0) ) > 1 ) trimmed1 = trimmed1[-tooextreme1 ]
tmp0[imeans] = mean( tmp0[ trimmed0 ], na.rm=TRUE )
tmp1[imeans] = mean( tmp1[ trimmed1 ], na.rm=TRUE )
O$bottom0 = NA
O$bottom1 = NA
O$bottom0.sd = NA
O$bottom1.sd = NA
O$bottom0.n = NA
O$bottom1.n = NA
O$bottom.diff = NA
O$bottom.diff.sd = NA
if ( any (is.na( c( tmp0[ standard ], tmp1[ standard ]) ) ) ) {
# no manual standard .. use mean as the standard
O$bottom0 = tmp0[imeans]
O$bottom1 = tmp1[imeans]
O$bottom0.sd = sd( ( tmp0[ trimmed0]), na.rm=TRUE ) # in secconds
O$bottom1.sd = sd( ( tmp1[ trimmed1]), na.rm=TRUE )
O$bottom0.n = length( which( is.finite( tmp0[ trimmed0] )) )
O$bottom1.n = length( which( is.finite( tmp1[ trimmed1] )) )
O$bottom.diff = difftime( O$bottom1, O$bottom0, units="secs" )
O$bottom.diff.sd = sqrt(O$bottom0.sd ^2 + O$bottom0.sd^2) # sec
} else {
# over-ride all data and use the manually determined results
O$bottom0 = tmp0[ standard ]
O$bottom1 = tmp1[ standard ]
O$bottom0.sd = NA
O$bottom1.sd = NA
O$bottom0.n = 1
O$bottom1.n = 1
O$bottom.diff = difftime( O$bottom1, O$bottom0, units="secs" )
}
tmp = data.frame( start=tmp0)
tmp$end = tmp1
tmp$diff = difftime( tmp[,"end"], tmp[,"start"], units="secs" )
tmp$start.bias = difftime( tmp[,"start"], O$bottom0, units="secs" )
tmp$end.bias = difftime( tmp[,"end"], O$bottom1, units="secs" )
rownames(tmp) = bcmethods
O$summary = tmp
# finalised data which have been filtered
fin.all = which( x$timestamp >= O$bottom0 & x$timestamp <= O$bottom1 )
if (length( fin.all ) == 0 ) fin.all = min( which( O$good)) : max( which( O$good) )
fin.good = intersect( fin.all, which( O$good) )
fin0 = min( fin.all, na.rm=TRUE)
fin1 = max( fin.all, na.rm=TRUE)
O$depth.mean = mean( x$depth[ fin.good ], na.rm=TRUE )
O$depth.sd = sd( x$depth[ fin.good ], na.rm=TRUE)
O$depth.n = length( fin.good )
O$depth.n.bad = length( fin.all) - O$depth.n
O$depth.smoothed.mean = mean( x$depth.smoothed[ fin0:fin1 ], na.rm=TRUE )
O$depth.smoothed.sd = sd( x$depth.smoothed[ fin0:fin1 ], na.rm=TRUE )
O$good = O$good
O$depth.filtered = fin.good
O$depth.smoothed = x$depth.smoothed
O$ts = x$ts
O$timestamp = x$timestamp
O$signal2noise = O$depth.n / length( fin.all ) # not really signal to noise but rather % informations
O$bottom.contact.indices = fin.all
O$bottom.contact = rep( FALSE, nx )
O$bottom.contact[ fin.all ] = TRUE
# estimate surface area if possible using wingspread &/or doorspread
O$surface.area = NA
if ( "wingspread" %in% names(x) | "doorspread" %in% names(x) ) {
sa = try( surfacearea.estimate( bcp=bcp, O=O ), silent=TRUE )
if ( ! "try-error" %in% class( sa ) ) O$surface.area = sa
}
# for minilog and seabird data .. we have temperature estimates to make ..
tmean= NA
tmeansd = NA
if ("temperature" %in% names(x) ) {
tmean = mean( x$temperature[fin.all], na.rm=TRUE )
tmeansd = sd( x$temperature[fin.all], na.rm=TRUE )
}
O$res = data.frame( cbind(z=O$depth.mean, t=tmean, zsd=O$depth.sd, tsd=tmeansd,
n=O$depth.n, t0=O$bottom0, t1=O$bottom1, dt=O$bottom.diff ) ) # this is really for the snow crab system
# time format / zone gets reset ..
O$res$t0 = as.POSIXct( O$res$t0, origin=lubridate::origin, tz="UTC" )
O$res$t1 = as.POSIXct( O$res$t1, origin=lubridate::origin, tz="UTC" )
O$res$dt = O$res$t1 - O$res$t0
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "yellow"
points( depth~ts, x[ O$maxdepth.method.indices, ], pch=20, col=mcol, cex=0.2)
}
print( O$summary)
return( O )
}
jae0/netmensuration documentation built on Jan. 11, 2024, 3:35 a.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 bottom.contact
Documented in bottom.contact
#' @title bottom.contact
#' @description Calculates indicies of touchdown and liftoff
#' @import tcltk geosphere lubridate
#' @return dataframe with various results
#' @export
bottom.contact = function( x, bcp, debugrun=FALSE ) {
require(tcltk)
# all timestamps must be in Posix/UTC
# range comparisons seem to fail when they are not
tzone = "UTC"
if (debugrun) {
debug.plot = TRUE
browser()
}
debug.plot = FALSE
n.min.required = 30
nx = nrow(x)
O = list() # output list
O$id = bcp$id
O$error.flag = NA
O$good = rep(TRUE, nx) # rows that will contain data that passes each step of data quality checks
if (length (which( (!is.na(x$depth)))) < n.min.required ) return( NULL )
##--------------------------------
# sort in case time is not in sequence
# timestamps have frequencies higher than 1 sec .. duplciates are created and this can pose a problem
x = x[order( x$timestamp ) ,]
x$ts = as.numeric( difftime( x$timestamp, min(x$timestamp), units="secs" ) )
O$plotdata = x # save incoming data after creation of time stamps and ordering
## PRE-FILTER 1
if ( exists( "double.depth.sensors", bcp) ) {
# two depth sensors were used simultaneously but they are not calibrated!
# remarkably hard to filter this out with any reliability at a later stage ..
om = modes( x$depth )
oaoi = range( which( x$depth <= om$ub2 & x$depth >= om$lb2))
oj = oaoi[1]:oaoi[2]
oi = interpolate.xy.robust( x[oj, c("ts", "depth")], method="loess", trim=0.05, probs=c(0.05, 0.95), target.r2=0.5 )
oo = which( (x$depth[oj] - oi) < 0)
x$depth[oj[oo]] = NA
if (any( is.na( x$depth))) x$depth = interpolate.xy.robust( x[, c("ts", "depth")], method="simple.linear" )
}
## PRE-FILTER 2
if (all(!is.na(bcp$time.gate))) {
# simple time-based gating if a time range is provided..
bcp$trange.max = (bcp$tdif.max+5)
O$good = bottom.contact.gating.time ( Zt=x$timestamp, good=O$good, bcp=bcp )
x$depth[ !O$good ] = NA
}
## PRE-FILTER 3
# simple depth-based gating
if( !any(x$depth > bcp$depth.min,na.rm=T)) return( NULL )
O$good = bottom.contact.gating.depth ( Z=x$depth, good=O$good, bcp=bcp)
x$depth[ !O$good ] = NA
## PRE-FILTER 4
# time and depth-based gating
mm = modes( x$depth ) # naive first estimate of location of most data depths
if (mm$lb == mm$ub || is.na(mm$sd)) return(NULL) # there is no depth variation in the data ... likely a bad data series
mm.i = which( x$depth > (mm$lb2+bcp$depth.range[1]) & x$depth < (mm$ub2 + bcp$depth.range[2]) )
O$good[ setdiff(1:nx, mm.i)] = FALSE
O$good[mm.i] = TRUE
## RANGE CHECKS
x.timerange = range( x$timestamp[O$good], na.rm=TRUE )
x.dt = difftime( x.timerange[2], x.timerange[1], units="mins" )
if ( x.dt < (bcp$tdif.min ) ) {
O$error.flag = "Time track is too short?"
return(O)
}
if ( x.dt > (bcp$tdif.max+10) ) { # +10 = 5 min on each side
# data vector is too long ... truncate
mm.r = range( mm.i )
mm.dt = difftime( x$timestamp[mm.r[2]], x$timestamp[mm.r[1]], units="mins" )
if ( mm.dt > (bcp$tdif.max+10) ) {
O$error.flag = "Time track is too long?"
}
}
if ( sd(x$depth, na.rm=TRUE) < bcp$eps.depth ) {
O$error.flag = "not enough variability in data?"
return(O)
}
## inSANITY CHECKS
# sometimes multiple tows exist in one track ...
# over-smooth depth to capture strange tows
x$dsm = interpolate.xy.robust( x[, c("ts", "depth")], method="sequential.linear", trim=0.05 )
x$dsm = interpolate.xy.robust( x[, c("ts", "dsm")], method="local.variance", trim=0.05 )
x$dsm = interpolate.xy.robust( x[, c("ts", "dsm")], method="sequential.linear", trim=0.05 )
x$dsm = interpolate.xy.robust( x[, c("ts", "dsm")], method="loess", trim=0.05 )
zz = rep( 0, nx )
zz[ which( x$dsm < ( mm$mode / 3 )) ] = 1 # flag shallow areas
inc.depth = abs( diff( x$dsm ) )
# also capture strong noise - very obviously wrong data
rapid.depth.changes = which( inc.depth > bcp$maxdepthchange )
if ( length( rapid.depth.changes ) > 0 ) {
zz[ rapid.depth.changes ] = 1
zz[ rapid.depth.changes-1 ] = 1 # include adjecent points to remove
zz[ rapid.depth.changes+1 ] = 1
}
dzz = diff(zz)
bnds = c(1, which( dzz != 0 ), nx )
if ( length (bnds) > 2 ) {
# i.e. , contaminated by noise or multiple tows
segs = diff(bnds) # length of each segment
longest = which.max( segs )
gg = bnds[longest]:bnds[(longest+1)]
bad = setdiff( 1:nx, gg)
O$good[bad] = FALSE
}
# fix tails
mm = modes( x$depth[ O$good ] ) # recompute after gating
mm.bot = which ( x$depth >= mm$lb2 & x$depth <= mm$ub2 ) # first est of bottom
if(length(mm.bot) == 0) return(NULL) # not enough variability in data
mm.med = floor( median( mm.bot, na.rm=TRUE ) )
# fix tails when the tails are not a single smooth tail:
llim = mm$mode + bcp$depth.range[1]
for ( ll in mm.med:1 ) {
if ( is.finite( x$depth[ll] ) ) {
if ( x$depth[ll] < llim ) break()
}
}
for ( uu in mm.med:nrow(x)) {
if ( is.finite( x$depth[uu] ) ) {
if ( x$depth[uu] < llim ) break()
}
}
todrop = c( 1:ll, uu:nrow(x) )
x$depth[ todrop ] = NA
O$good[ todrop] = FALSE
if ( sd(x$depth, na.rm=TRUE) < bcp$eps.depth ) return(NULL) # not enough variability in data
## ------------------------------
## MAIN NOISE/INTERPOLATION FILTER
# Some filtering of noise from data and further focus upon area of interest based upon time and depth if possible
res = NULL
res = try( bottom.contact.filter.noise ( x=x, good=O$good, bcp=bcp, debug=debugrun ), silent =TRUE )
if ( "try-error" %in% class(res) || is.null(res$depth.smoothed) ) {
x$depth = jitter( x$depth )
res = try( bottom.contact.filter.noise ( x=x, good=O$good, bcp=bcp, debug=debugrun ), silent =TRUE )
}
if ( !"try-error" %in% class(res) ) {
if (!is.null(res$depth.smoothed)) {
if ( cor( x$depth, res$depth.smoothed, use="pairwise.complete.obs") > 0.999 || is.na(cor( x$depth, res$depth.smoothed, use="pairwise.complete.obs") )) {
bcp$noisefilter.target.r2 = bcp$noisefilter.target.r2 - 0.1
x$depth = jitter( x$depth )
res = try( bottom.contact.filter.noise ( x=x, good=O$good, bcp=bcp, debug=debugrun ), silent =TRUE )
}}}
x$depth.smoothed = x$depth
if ( ! "try-error" %in% class( res) ) {
# retain good/bad information only for the interval where there is data
# .. this effectively does the filtering in the area of interest only and retains the tails
# for further influence later
ig = range( which( res$good ))
igg = ig[1]:ig[2]
O$good[igg] = res$good[igg]
x$depth[ !O$good ] = NA
if (!is.null(res$depth.smoothed)) x$depth.smoothed= res$depth.smoothed
}
if(sum(x$depth-min(x$depth,na.rm=T),na.rm=T)==0) return( NULL )
if(sum(O$good)==0) return( NULL )
if (sd(x$depth, na.rm=TRUE) < bcp$eps.depth ) return(NULL) # not enough variability in data
if (sd(x$depth.smoothed, na.rm=TRUE) < bcp$eps.depth ) return(NULL) # not enough variability in data
x.timerange = range( x$timestamp[O$good], na.rm=TRUE )
x.dt = difftime( x.timerange[2], x.timerange[1], units="mins" )
if ( x.dt < ( bcp$tdif.min ) ) {
O$error.flag = "Not enough data?"
return(O)
}
if ( x.dt > (bcp$tdif.max + 10 ) ) {
O$error.flag = "Too much data?"
return(O)
}
## FINAL DATA GATING
# variance gating attempt after the data has been cleaned as much as possible
O$variance.method0 = NA
O$variance.method1 = NA
O$variance.method.indices = NA
res = NULL
res = try( bottom.contact.gating.variance ( x, O$good, bcp ), silent =TRUE )
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1 ) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$variance.method0 = res$bc0
O$variance.method1 = res$bc1
O$variance.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 )
bad = which( x$timestamp < res$bc0 | x$timestamp > res$bc1 )
if (length( bad) > 0) O$good[ bad ] = FALSE
x$depth[ !O$good ] = NA
} }
}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "gray"
points( depth~ts, x[ O$variance.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## NOTE::: From this point on, O$good is now complete
## -- it contains indices of time and depth-based gating as well as the variance based gating
# finalize selection of area of interest (based upon gating, above)
aoi.range = range( which( O$good ) )
aoi.mid = floor( mean( aoi.range ) ) # approximate midpoint
aoi.min = aoi.range[1]
aoi.max = aoi.range[2]
O$aoi = aoi.min:aoi.max # stored for use in the linear method where we need to recover the left-most index
##--------------------------------
# Modal method: gating by looking for modal distribution and estimating sd of the modal group in the data
# first by removing small densities ( 1/(length(i)/nb) ) and by varying the number of breaks in the histogram
# until a target number of breaks, nbins with valid data are found
O$modal.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$modal.method1 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$modal.method.indices = NA
sm0 = x[ O$aoi, c("depth.smoothed", "timestamp", "ts" ) ] # send filtered data ... continuity not important .. order is important
res = NULL
res = try( bottom.contact.modal( sm=sm0, bcp ), silent=TRUE )
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1 ) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$modal.method0 = res$bc0 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$modal.method1 = res$bc1 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$modal.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 ) # x correct
} }
}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "green"
points( depth~ts, x[ O$modal.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## ----------------------------
## Smooth method: using smoothed data (slopes are too unstable with raw data),
## compute first derivatives to determine when the slopes inflect
O$smooth.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$smooth.method1 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$smooth.method.indices = NA
sm0 = x[ O$aoi, c("depth.smoothed", "timestamp", "ts")] # Send all data within the aoi --- check this .. order is important
res = NULL
res = try(
bottom.contact.smooth( sm=sm0, bcp=bcp ) , silent =TRUE)
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$smooth.method0 = res$bc0 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$smooth.method1 = res$bc1 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$smooth.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 ) # x correct
} }
}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "blue"
points( depth~ts, x[ O$smooth.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## ----------------------------
## maxdepth method: looking for the max depth near the areas of interest (left, right)
O$maxdepth.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$maxdepth.method1 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$maxdepth.method.indices = NA
sm0 = x[, c("depth.smoothed", "timestamp", "ts")] # Send all data within the aoi --- check this .. order is important
sm0$depth[ !O$good ] = NA
sm0 = sm0[ O$aoi, ]
bcmethods=c( "smooth.method", "modal.method" )
res = NULL
res = try( bottom.contact.maxdepth( sm=sm0, O=O, bcmethods=bcmethods, bcp=bcp ) , silent=TRUE )
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if(length(res$bc1)>0 && length(res$bc0)>0){
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1 ) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$maxdepth.method0 = res$bc0 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$maxdepth.method1 = res$bc1 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$maxdepth.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 ) # x correct
}}
}}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "yellow"
points( depth~ts, x[ O$maxdepth.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## ---------------------------
## Linear method: looking at the intersection of three lines (up, bot and down)
O$linear.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$linear.method1 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$linear.method.indices = NA
sm0 = x[, c("depth.smoothed", "timestamp", "ts")] # Send all data within the aoi --- check this .. order is important
sm0$depth[ !O$good ] = NA
sm0 = sm0[ O$aoi, ]
bcmethods=c( "smooth.method", "modal.method", "linear.method" )
res = NULL
res = try( bottom.contact.linear( sm=sm0, O=O, bcmethods=bcmethods, bcp=bcp ) , silent=TRUE )
if ( ! "try-error" %in% class( res) ) {
if ("bc0" %in% names(res)) {
if ( all(is.finite( c(res$bc0, res$bc1 )) ) ) {
DT = abs( as.numeric( difftime( res$bc0, res$bc1, units="mins" ) ) )
if ( length(DT) == 1) {
if ( is.finite(DT) && DT > bcp$tdif.min & DT < bcp$tdif.max ) {
O$linear.method0 = res$bc0 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$linear.method1 = res$bc1 #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$linear.method.indices = which( x$timestamp >= res$bc0 & x$timestamp <= res$bc1 ) # x correct
} }
}}
}
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "red"
points( depth~ts, x[ O$linear.method.indices, ], pch=20, col=mcol, cex=0.2)
}
## ---------------------------
O$manual.method0 = NA #### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$manual.method1 = NA ### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
if ( bcp$user.interaction ) {
satisified = F
while(!satisified){
#open plot window, user can modify the numbers to fit their screen
x11(width = 25, height = 15)
#Define Margins for multiple y-axis
par(mar=c(5, 8, 4, 4) + 0.1)
grange = min(which(O$good)):max(which(O$good))
if("wingspread" %in% names(O$plotdata)){
#Plot the door spread
ylim=c(min(O$plotdata$wingspread[grange], na.rm = TRUE)-1,max(O$plotdata$wingspread[grange], na.rm = TRUE)+1)
if(is.finite(ylim)){
plot(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$wingspread[grange]))], O$plotdata$wingspread[grange][which(!is.na(O$plotdata$wingspread[grange]))], axes=F, ylim=ylim, xlab="", ylab="",type="p",col="#0000FF1A", main="",xlim=c(min(O$plotdata$timestamp[grange]), max(O$plotdata$timestamp[grange])))
if(length(unique(na.omit(O$plotdata$wingspread[grange]))) > 5 ){
smoothingSpline = smooth.spline(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$wingspread[grange]))], O$plotdata$wingspread[grange][which(!is.na(O$plotdata$wingspread[grange]))], spar=.5)
lines(smoothingSpline, col="blue")
}
abline( h = c(median(O$plotdata$wingspread[grange], na.rm = TRUE)), col = "blue", lty =2 )
axis(2,col="blue", col.lab = "blue", col.axis = "blue", lwd=1)
mtext(2,text="Wing Spread",col = "blue",line=2)
}
}
yl=c(min(O$plotdata$opening[grange], na.rm = TRUE)-1,max(O$plotdata$opening[grange], na.rm = TRUE)+1)
if(!is.numeric(yl[1]) | is.na(yl[1]) ) yl[1] = 0
if(!is.numeric(yl[2]) | is.na(yl[2]) ) yl[2] = 10
#Plot the opening
if("opening" %in% names(O$plotdata)){
if(length(which(!is.na(O$plotdata$opening[grange]))) > 10){
par(new=T)
plot(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$opening[grange]))], O$plotdata$opening[grange][which(!is.na(O$plotdata$opening[grange]))], axes=F, ylim=yl, xlab="", ylab="",type="p",col="#A52A2A1A", main="",xlim=c(min(O$plotdata$timestamp[grange]), max(O$plotdata$timestamp[grange])))
smoothingSpline = smooth.spline(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$opening[grange]))], O$plotdata$opening[grange][which(!is.na(O$plotdata$opening[grange]))], spar=.5)
lines(smoothingSpline, col="brown")
axis(2, ylim=yl,col = "brown",col.lab = "brown", col.axis = "brown",lwd=1,line=3.5)
mtext(2,text="Opening", col = "brown", line=5.5)
}
}
#Plot the depth (Clickable)
par(new=T)
plot(O$plotdata$timestamp[grange][which(!is.na(O$plotdata$depth[grange]))], O$plotdata$depth[grange][which(!is.na(O$plotdata$depth[grange]))], axes=F, ylim=rev(c(min(O$plotdata$depth[grange], na.rm = TRUE)-1,max(O$plotdata$depth[grange], na.rm = TRUE)+1)), xlab="", ylab="", type="l", main=O$id,xlim=c(min(O$plotdata$timestamp[grange]), max(O$plotdata$timestamp[grange])),lwd=1)
axis(4, ylim=rev(c(min(O$plotdata$depth[grange], na.rm = TRUE)-1,max(O$plotdata$depth[grange], na.rm = TRUE)+1)),lwd=1)
mtext(4,text="Depth", line = 2)
#Set up time axis plotting
abli = seq(O$plotdata$timestamp[grange][1],O$plotdata$timestamp[grange][length(O$plotdata$timestamp[grange])], "1 min")
#Draw the time axis
axis.POSIXct(1, at = abli, format = "%H:%M:%S", labels = TRUE)
ablit = unique(as.character(lubridate::date(abli)))
mtext(text = ablit, side=1, col="black", line=2)
abline( v = abli, col = "lightgrey")
#Visualize neutral haulback
#Find positions of going backward toward origin. Nuetral haulback
disfromorigin = rep(0, 1, length(O$plotdata$latitude[grange]))
for(k in 1:(length(O$plotdata$latitude[grange])-1)){
p1 = c(O$plotdata$longitude[grange][1],O$plotdata$latitude[grange][1])
p2 = c(O$plotdata$longitude[grange][k+1],O$plotdata$latitude[grange][k+1])
if(!is.null(p1) & !is.null(p2)){
if(!is.na(p1) & !is.na(p2))
{
disfromorigin[k+1] = geosphere::distHaversine(p1 , p2, r=6378137)
}
}
}
neu = rep(0, 1, length(O$plotdata$latitude[grange]))
for(k in 1:(length(disfromorigin)-1)){
neu[k] = disfromorigin[k+1]-disfromorigin[k]
if(k > 1){
if(neu[k-1]<0 & neu[k]>=0) neu[k-1] = 0
}
}
neuhaul = F
indn = which(neu < 0)
if(length(indn) > 10){
print(paste("Neutral Haulback! ", O$id))
neuhaul = T
}
if(length(indn)>0){
abline(v = O$plotdata$timestamp[grange][indn], col = "red")
}
id = identify(O$plotdata$timestamp[grange], O$plotdata$depth[grange], labels = c(as.character(O$plotdata$timestamp[grange])),n = 1, pos = TRUE)
start = O$plotdata$timestamp[grange][id$ind]
points(O$plotdata$timestamp[grange][id$ind], O$plotdata$depth[grange][id$ind], col = "darkgreen", pch = 19)
#Code that allow identifying clicks
id2 = identify(O$plotdata$timestamp[grange], O$plotdata$depth[grange], labels = c(as.character(O$plotdata$timestamp[grange])),n = 1, pos = TRUE)
end = O$plotdata$timestamp[grange][id2$ind]
points(O$plotdata$timestamp[grange][id2$ind], O$plotdata$depth[grange][id2$ind], col = "darkgreen", pch = 19)
##Code taken from bottom.contact.plot.r
legendtext = NULL
legendcol = NULL
legendpch = NULL
if (all(is.finite(c( O$variance.method0, O$variance.method1) ) ) ) {
mcol = "pink"
# points( depth~ts, x[ O$variance.method.indices, ], pch=20, col=mcol, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$variance.method.indices)], col=mcol, lty="solid", lwd=1.2)
abline (v=O$plotdata$timestamp[max(O$variance.method.indices)], col=mcol, lty="solid", lwd=1.2)
DT = as.numeric( difftime( O$variance.method1, O$variance.method0, units="mins" ) )
legendtext = c( legendtext, paste( "variance gate: ", round( DT, 2), "" ) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20 )
}
if (all(is.finite( c(O$modal.method0, O$modal.method1 ) ) ) ) {
mcol = "red" # colour for plotting
# points( depth~ts, x[O$modal.method.indices,], col=mcol, pch=20, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$modal.method.indices)], col=mcol, lty="dashed")
abline (v=O$plotdata$timestamp[max(O$modal.method.indices)], col=mcol, lty="dashed")
DT = as.numeric( difftime( O$modal.method1, O$modal.method0, units="mins" ) )
legendtext = c( legendtext, paste( "modal: ", round( DT, 2) ) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20)
}
if ( all(is.finite( c( O$smooth.method0, O$smooth.method1) ) ) ) {
mcol = "blue"
# points( depth~ts, x[O$smooth.method.indices,], col=mcol, pch=20, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$smooth.method.indices)], col=mcol, lty="dotdash", lwd=1.5)
abline (v=O$plotdata$timestamp[max(O$smooth.method.indices)], col=mcol, lty="dotdash", lwd=1.5)
DT = as.numeric( difftime( O$smooth.method1, O$smooth.method0, units="mins" ) )
legendtext = c(legendtext, paste( "smooth: ", round(DT, 2)) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20)
}
if (all(is.finite( c(O$linear.method0, O$linear.method1) )) ) {
mcol ="green"
# points( depth~ts, x[O$linear.method.indices,], col=mcol, pch=20, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$linear.method.indices)], col=mcol, lty="twodash")
abline (v=O$plotdata$timestamp[max(O$linear.method.indices)], col=mcol, lty="twodash")
DT = as.numeric( difftime( O$linear.method1, O$linear.method0, units="mins" ) )
legendtext = c( legendtext, paste( "linear: ", round( DT, 2) ) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20)
}
if (all(is.finite( c( O$maxdepth.method0, O$maxdepth.method1) )) ) {
mcol ="orange"
# points( depth~ts, x[O$linear.method.indices,], col=mcol, pch=20, cex=0.2)
abline (v=O$plotdata$timestamp[min(O$maxdepth.method.indices)], col=mcol, lty="solid")
abline (v=O$plotdata$timestamp[max(O$maxdepth.method.indices)], col=mcol, lty="solid")
DT = as.numeric( difftime( O$maxdepth.method1, O$maxdepth.method0, units="mins" ) )
legendtext = c( legendtext, paste( "maxdepth: ", round( DT, 2) ) )
legendcol = c( legendcol, mcol)
legendpch =c( legendpch, 20)
}
if ( !( is.null( legendtext))) legend( "top", legend=legendtext, col=legendcol, pch=legendpch, bg = "white" )
#END code taken from bottom.contact
#User interface
goodans = FALSE
res <- tkmessageBox(title = "ManualTD",
message = "Are you happy with these values?", icon = "info", type = "yesno")
if(grepl("yes", res)){
satisified = TRUE
goodans = TRUE
O$manual.method0 = start
O$manual.method1 = end
O$manual.method.indices = which( O$plotdata$timestamp >= O$manual.method0 & O$plotdata$timestamp <= O$manual.method1 )
# result = rbind(result, row)
# if(!file.exists(file)){
# write.table(row, file = file, sep = ",", row.names = F, quote = F)
# }else{
# write.table(row, file = file, sep = ",", append = T, col.names = F, row.names = F, quote = F )
# }
}
if(grepl("no", res)){
satisified = FALSE
goodans = TRUE
}
if(!goodans)print("Did not receive valid command, please try again.")
dev.off()
}#END plotting loop while not satisfied
if(is.na(O$manual.method0) || is.na(O$manual.method1)){
bcp$user.interaction = FALSE
}
else{
if(is.null(bcp$station)) station = unlist(strsplit(bcp$id, "\\."))[4]
else station = bcp$station
bcp$id = bcp$trip
mf = NULL
if(!is.null(bcp$from.manual.archive)) mf = file.path(bcp$from.manual.archive, "clicktouchdown_all.csv")
if(!is.null(bcp$from.manual.file)) mf = bcp$from.manual.file
if(!is.null(mf)){
manualclick = NULL
if(file.exists(mf)){
manualclick = read.csv(mf, as.is=TRUE)
if(bcp$datasource == "lobster"){
sta.ind = which(manualclick$station == station & manualclick$trip == bcp$trip)
}
else{
sta.ind = which(manualclick$station == station & manualclick$year == bcp$YR)
}
if(length(sta.ind) == 0){
print(bcp$trip)
manualclick = rbind(manualclick, data.frame(station = station, start = unlist(as.character(O$manual.method0)), end = unlist(as.character(O$manual.method1)), depth = mean( x$depth, na.rm=TRUE ), year = bcp$YR, trip = bcp$id))
}
else{
manualclick$station[sta.ind] = station
manualclick$start[sta.ind] = as.character(O$manual.method0)
manualclick$end[sta.ind] = as.character(O$manual.method1)
manualclick$depth[sta.ind] = mean( x$depth, na.rm=TRUE )
manualclick$year[sta.ind] = bcp$YR
manualclick$trip[sta.ind] = bcp$id
}
}
else{
manualclick = data.frame(station = station, start = as.character(O$manual.method0), end = as.character(O$manual.method1), depth = mean( x$depth, na.rm=TRUE ), year = bcp$YR, trip = bcp$id)
}
write.csv(manualclick, mf, row.names = FALSE )
}
}
}
#BC- Added condition in case user interaction is desired, no need to do this step.
if (!is.null(bcp$from.manual.archive) && !bcp$user.interaction) {
print( "Loading values from previously generated .csv")
if(file.exists(file.path(bcp$from.manual.archive, "clicktouchdown_all.csv"))){
manualclick = read.csv(file.path(bcp$from.manual.archive, "clicktouchdown_all.csv"), as.is=TRUE)
station = unlist(strsplit(bcp$id, "\\."))[4]
sta.ind = which(manualclick$station == station & manualclick$year == bcp$YR)
if(length(sta.ind == 1)){
O$manual.method0 = lubridate::ymd_hms(manualclick$start[sta.ind], tz = "UTC")
O$manual.method1 = lubridate::ymd_hms(manualclick$end[sta.ind], tz = "UTC")
}
}
}
O$means0 = NA ### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
O$means1 = NA ### NOTE:: using the 'c' operator on posix strips out the timezone info! this must be retained
bcmethods = c("manual.method", "variance.method", "smooth.method", "modal.method", "maxdepth.method", "linear.method", "means" )
standard = which( bcmethods=="manual.method") # gold standard
direct = which( bcmethods %in% c("smooth.method", "modal.method", "linear.method", "maxdepth.method" ) )
imeans = which( bcmethods == "means" )
# must be careful as timestamp is being converted to text and will lose tzone ... need to reset to correct timezone:
bcm0 = paste(bcmethods, "0", sep="")
bcm1 = paste(bcmethods, "1", sep="")
# recovert to time zone of incoming data as time zone is lost with the transpose
tmp0 = lubridate::ymd_hms( t( as.data.frame(O[ bcm0 ]) ) ) # UTC
tmp1 = lubridate::ymd_hms( t( as.data.frame(O[ bcm1 ]) ) ) # UTC
bottom0.mean = mean(tmp0, na.rm=TRUE)
bottom1.mean = mean(tmp1, na.rm=TRUE)
bottom0.sd = max( 1, sd( tmp0, na.rm=TRUE ) ) # in seconds
bottom1.sd = max( 1, sd( tmp1, na.rm=TRUE ) )
dflag0 = rowSums(as.matrix(dist(tmp0[direct], upper=TRUE )), na.rm=TRUE) # which are most extreme
dflag1 = rowSums(as.matrix(dist(tmp1[direct], upper=TRUE )), na.rm=TRUE) # which are most extreme
tooextreme0 = which.max( dflag0 )
tooextreme1 = which.max( dflag1 )
trimmed0 = trimmed1 = direct
if (length( which(dflag0 > 0) ) > 1 ) trimmed0 = trimmed0[-tooextreme0 ]
if (length( which(dflag1 > 0) ) > 1 ) trimmed1 = trimmed1[-tooextreme1 ]
tmp0[imeans] = mean( tmp0[ trimmed0 ], na.rm=TRUE )
tmp1[imeans] = mean( tmp1[ trimmed1 ], na.rm=TRUE )
O$bottom0 = NA
O$bottom1 = NA
O$bottom0.sd = NA
O$bottom1.sd = NA
O$bottom0.n = NA
O$bottom1.n = NA
O$bottom.diff = NA
O$bottom.diff.sd = NA
if ( any (is.na( c( tmp0[ standard ], tmp1[ standard ]) ) ) ) {
# no manual standard .. use mean as the standard
O$bottom0 = tmp0[imeans]
O$bottom1 = tmp1[imeans]
O$bottom0.sd = sd( ( tmp0[ trimmed0]), na.rm=TRUE ) # in secconds
O$bottom1.sd = sd( ( tmp1[ trimmed1]), na.rm=TRUE )
O$bottom0.n = length( which( is.finite( tmp0[ trimmed0] )) )
O$bottom1.n = length( which( is.finite( tmp1[ trimmed1] )) )
O$bottom.diff = difftime( O$bottom1, O$bottom0, units="secs" )
O$bottom.diff.sd = sqrt(O$bottom0.sd ^2 + O$bottom0.sd^2) # sec
} else {
# over-ride all data and use the manually determined results
O$bottom0 = tmp0[ standard ]
O$bottom1 = tmp1[ standard ]
O$bottom0.sd = NA
O$bottom1.sd = NA
O$bottom0.n = 1
O$bottom1.n = 1
O$bottom.diff = difftime( O$bottom1, O$bottom0, units="secs" )
}
tmp = data.frame( start=tmp0)
tmp$end = tmp1
tmp$diff = difftime( tmp[,"end"], tmp[,"start"], units="secs" )
tmp$start.bias = difftime( tmp[,"start"], O$bottom0, units="secs" )
tmp$end.bias = difftime( tmp[,"end"], O$bottom1, units="secs" )
rownames(tmp) = bcmethods
O$summary = tmp
# finalised data which have been filtered
fin.all = which( x$timestamp >= O$bottom0 & x$timestamp <= O$bottom1 )
if (length( fin.all ) == 0 ) fin.all = min( which( O$good)) : max( which( O$good) )
fin.good = intersect( fin.all, which( O$good) )
fin0 = min( fin.all, na.rm=TRUE)
fin1 = max( fin.all, na.rm=TRUE)
O$depth.mean = mean( x$depth[ fin.good ], na.rm=TRUE )
O$depth.sd = sd( x$depth[ fin.good ], na.rm=TRUE)
O$depth.n = length( fin.good )
O$depth.n.bad = length( fin.all) - O$depth.n
O$depth.smoothed.mean = mean( x$depth.smoothed[ fin0:fin1 ], na.rm=TRUE )
O$depth.smoothed.sd = sd( x$depth.smoothed[ fin0:fin1 ], na.rm=TRUE )
O$good = O$good
O$depth.filtered = fin.good
O$depth.smoothed = x$depth.smoothed
O$ts = x$ts
O$timestamp = x$timestamp
O$signal2noise = O$depth.n / length( fin.all ) # not really signal to noise but rather % informations
O$bottom.contact.indices = fin.all
O$bottom.contact = rep( FALSE, nx )
O$bottom.contact[ fin.all ] = TRUE
# estimate surface area if possible using wingspread &/or doorspread
O$surface.area = NA
if ( "wingspread" %in% names(x) | "doorspread" %in% names(x) ) {
sa = try( surfacearea.estimate( bcp=bcp, O=O ), silent=TRUE )
if ( ! "try-error" %in% class( sa ) ) O$surface.area = sa
}
# for minilog and seabird data .. we have temperature estimates to make ..
tmean= NA
tmeansd = NA
if ("temperature" %in% names(x) ) {
tmean = mean( x$temperature[fin.all], na.rm=TRUE )
tmeansd = sd( x$temperature[fin.all], na.rm=TRUE )
}
O$res = data.frame( cbind(z=O$depth.mean, t=tmean, zsd=O$depth.sd, tsd=tmeansd,
n=O$depth.n, t0=O$bottom0, t1=O$bottom1, dt=O$bottom.diff ) ) # this is really for the snow crab system
# time format / zone gets reset ..
O$res$t0 = as.POSIXct( O$res$t0, origin=lubridate::origin, tz="UTC" )
O$res$t1 = as.POSIXct( O$res$t1, origin=lubridate::origin, tz="UTC" )
O$res$dt = O$res$t1 - O$res$t0
if(debug.plot) {
trange = range( x$ts[O$good], na.rm=TRUE )
drange = c( quantile( x$depth, c(0.05, 0.975), na.rm=TRUE) , median( x$depth, na.rm=TRUE ) * 1.05 )
#BC - Plots fail in RStudio graphics device, add clause
dev.new(noRStudioGD = TRUE)
plot(depth~ts, x, ylim=c(drange[2],drange[1]), xlim=c(trange[1],trange[2]), pch=20, cex=0.1, col="gray" )
mcol = "yellow"
points( depth~ts, x[ O$maxdepth.method.indices, ], pch=20, col=mcol, cex=0.2)
}
print( O$summary)
return( O )
}
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.