R/net.configuration.r
In jae0/bio.snowcrab: Snow crab assessment suite using aegis* and associated projects and data streams.
Defines functions
net.configuration
net.configuration = function( N, t0=NULL, t1=NULL, set_timestamp=NULL, yr=NULL, plotdata=TRUE ) {
# N is netmind data
# t0 is current best estimate of start and end time
# set_timestamp is timestamp from set-log .. alternate if there is no other useful time marker
# note:: NETMIND data stream calls the wingspread as "DOORSPREAD".
# It is not. It is wingspread.
plotdir = project.datadirectory("bio.snowcrab", "data", "netmind", "figures" )
# create default output should the following fail
out = data.frame( slon=NA, slat=NA, distance=NA, spread=NA, spread_sd=NA,
surfacearea=NA, vel=NA, vel_sd=NA, netmind_n=NA, t0=NA, t1=NA, dt=NA, yr=NA,
temperature.n=NA, temperature_sd.n=NA )
n.req = 10
if (!is.null(t0)) {
if (length(t0)>1) {t0 = NULL}
if (is.na(t0)) t0 = NULL
}
if (!is.null(t1)) {
if (is.na(t1)) t1 = NULL
}
#changed this switch from depgth filed to lat as if there is not depth info can still run script
if ( length( which( is.finite( N$lat))) < n.req ) print(N[1,])
problem = F
# time checks
if ( is.null(t0) ){
if ( !is.null(set_timestamp) ) {
t0 = set_timestamp # no data in t0 ,, use set_timestamp as alternate
set_timestamp =NULL # remove to cause no more effects
}
}
#bad.list = c('netmind.S26092014.9.541.17.48.304',
# 'netmind.S20092007.8.333.17.27.241' )
bad.list = NULL
bad.list = unique( c(bad.list, p$netmensuration.problems) )
if ( N$netmind_uid[1] %in% bad.list) {
return(out)
}
if ( any( is.null( t1 ) || is.null(t0) ) ) {
# try to determine from netmind data if minilog/seadbird data methods have failed. .. not effective due to noise/and small data stream
print(N$netmind_uid[1])
M = N[, c("timestamp", "depth") ]
oo = which( is.finite( M$depth ))
if ( length(oo) < n.req ) return(out)
if(!is.null(set_timestamp)) settimestamp=set_timestamp
if(is.null(set_timestamp)) settimestamp=t0
time.gate = list( t0=settimestamp - dminutes(6), t1=settimestamp + dminutes(12) )
bcp = list(id=N$netmind_uid[1], nr=nrow(M), YR=yr, trip = unlist(strsplit(N$netmind_uid[1], "\\."))[2], tdif.min=3, tdif.max=11, time.gate=time.gate,
depth.min=20, depth.range=c(-25,15), eps.depth=3,
smooth.windowsize=5, modal.windowsize=5, datasource = "snowcrab",
noisefilter.trim=0.05, noisefilter.target.r2=0.8, noisefilter.quants=c(0.05, 0.95) )
if(yr<2007)bcp$from.manual.archive=FALSE # manual touchdown only done since 2007
bcp = bottom.contact.parameters( bcp ) # add other default parameters
#BC: Determine if this station was done yet, if not then we want user interaction.
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) bcp$user.interaction = FALSE
else bcp$user.interaction = TRUE
}
dp.out.range = which(M$depth < .1 | M$depth > 360)
if(length(dp.out.range)>0){
M$depth[dp.out.range] = NA
}
M$wingspread = M$doorspread
bc = NULL
bc = bottom.contact( x=M, bcp=bcp )
if ( is.null(bc) || (!is.null( bc$res) && ( ( !is.finite(bc$res$t0 ) || !is.finite(bc$res$t1 ) ) ) )) {
bcp$noisefilter.target.r2=0.8
bc = bottom.contact( x=M, bcp=bcp )
}
if ( is.null(bc) || (!is.null( bc$res) && ( ( !is.finite(bc$res$t0 ) || !is.finite(bc$res$t1 ) ) ) )) {
# try once more with random settings
bcp$noisefilter.inla.h = 0.05
bcp$noisefilter.target.r2= 0.75
bc = bottom.contact( x=M, bcp=bcp )
}
if ( is.null(bc) || (!is.null( bc$res) && ( ( !is.finite(bc$res$t0 ) || !is.finite(bc$res$t1 ) ) ) )) {
# try once more with more random settings .. noise is high in netminds data
bcp$eps.depth = 5
M$depth = jitter( M$depth, amount = bcp$eps.depth/10 )
bcp$noisefilter.inla.h = 0.1
bc = bottom.contact( x=M, bcp=bcp )
}
if ( !is.null(bc) ) {
if (plotdata) {
bottom.contact.plot( bc )
plotfn = file.path( plotdir, paste(bcp$id, "pdf", sep="." ) )
print (plotfn)
dev.flush()
dev.copy2pdf( file=plotfn )
graphics.off()
}
}
if (is.null(t0) & !is.null(bc$bottom0) ) t0 = bc$bottom0
if (is.null(t1) & !is.null(bc$bottom1) ) t1 = bc$bottom1
N = N[ bc$bottom.contact , ]
}
if (all(is.na(t0))) t0=NA
if (all(is.na(t1))) t1=NA
if ( is.null(t1) || !is.finite(t1) ) {
t1_tmp = NA # rather than guessing, flag and then fill later
} else {
t1_tmp = t1
}
# if we are here, it is because a reasonable amount of data is present ..
# do some more checks and get a first estimate of some parameters in case other errors/limits are found
out$slon=N$lon[1]
out$slat=N$lat[1]
out$spread=mean( N$doorspread, na.rm=T ) / 1000 # though called doorspread it is actually wingspread
out$spread_sd=sd( N$doorspread, na.rm=T ) /1000
# first set the impossible door spreads to NA and then interpolate based upon a 5 and 95 % quantiles
# Netmind has very noisy data
hl = c( 3, 16 ) # m
ihl = which( N$doorspread < hl[1] | N$doorspread > hl[2] )
if (length (ihl)>0 ) N$doorspread[ihl] = NA
quantiles.to.trim = c(0.05, 0.95)
qnt = quantile( N$doorspread[N$doorspread>0], quantiles.to.trim, na.rm=T)
iqnt = which( N$doorspread < qnt[1] | N$doorspread > qnt[2] )
N$doorspread[ iqnt ] = NA
gooddoor = which( is.finite(N$doorspread) )
baddoor = which( ! is.finite(N$doorspread) )
if (length(baddoor)>0) N$doorspread[ baddoor ] = NA
if ( length( gooddoor) < n.req ) {
#problem =T turned this off December 23, 2013 as if we have decent gps data now we can calc distance
if ( length( gooddoor) >0 ) {
out$netmind_n <- length(gooddoor)
out$spread = mean( N$doorspread[ gooddoor ], na.rm=T ) / 1000
out$spread_sd = sd( N$doorspread[ gooddoor ], na.rm=T ) / 1000
}
}
if ( length(t0)>1) t0 = NA
out$t0 = t0
out$t1 = t1
out$dt = difftime( as.POSIXct(t1), as.POSIXct(t0) ) # minutes
if(!is.na(out$t0)) out$yr = lubridate::year( out$t0)
if(is.na(out$t0)) out$yr = lubridate::year(N$timestamp[1])
itime = which( N$timestamp >= t0 & N$timestamp <= t1 )
if ( length( itime) < n.req ) problem = T
if (problem) {
## Need to check on these at some point
out$t0 = t0
out$t1 = NA
out$dt = NA
return(out)
}
# eOR checks
# Higher GPS data accuracy in 2019 resulted in high approximation of distance travelled. This was due to frequent occasions of backwards
# motion (assuming from boat roll and pitch) which was added to the total distance. When looking at previous years this also seem to have
# but not to the same extent however removing the backwards motion would result in better data accuracy from year using e-sonar 2014+.
# The solution is to trim the data by 3/4, or, with the current recording rate of every 4 seconds instead of every 1 second.
# The above works when the metrics are repeated but not when we actually only record the values when we get a sensor hit. This does not
# work for the 2021 data. Trying to smooth latitude and longitude so that we dont have data loss.
# This again did not work as paths were way out in some cases. Resolved by applying a windowed average smoothing
# Keeping code for st_simplify in case
itime.smoothing = which( N$timestamp >= (t0 - lubridate::seconds(10)) & N$timestamp <= (t1 + lubridate::seconds(10) ))
rolling.ave = 6
Ntemp = N[ itime.smoothing,]
Ntemp$nlat = NA
Ntemp$nlong = NA
for(i in 1:nrow(Ntemp)){
inds = (i-rolling.ave):(i+rolling.ave)
inds = inds[which(inds>0)]
inds = inds[which(inds<=nrow(Ntemp))]
Ntemp$nlong[i] = mean(Ntemp$lon[inds])
Ntemp$nlat[i] = mean(Ntemp$lat[inds])
}
#plot(Ntemp$lon,Ntemp$lat, type = "l", col = "#000000E6", lwd = 3)
#lines(Ntemp$nlon,Ntemp$nlat, col = 'red', lwd = 2)
itime = which( Ntemp$timestamp >= t0 & Ntemp$timestamp <= t1 )
N = Ntemp[ itime, ]
if(yr>=2014){
# require(sf)
# tolerence = 5
# N$new.lat = NA
# N$new.lon = NA
# xy = st_cast( st_multipoint(cbind(N$lon,N$lat), dim = "XY"), 'LINESTRING' )
# cxyll = st_sfc(xy, crs = "+init=epsg:4326" )
# c_xy = st_transform(cxyll, st_crs("+proj=utm +zone=20 +ellps=GRS80 +datum=NAD83 +units=m" ))
# N$utm.y = c_xy[[1]][,2]
# N$utm.x = c_xy[[1]][,1]
# xy_smoothed = st_simplify( c_xy, dTolerance=tolerence)
# tol.count = tolerence+1
# while(any(grepl("MULTILINESTRING", class(xy_smoothed)))){
# tol.count = tol.count+1
# print(paste("Moving tolerence to:", tol.count))
# xy_smoothed = st_simplify( c_xy, dTolerance=tol.count)
# if(tol.count > 100) break()
# }
#
#
#
# nodes = which(paste(N$utm.y, N$utm.x) %in% paste(xy_smoothed[[1]][,2],xy_smoothed[[1]][,1]))
# nodes = nodes[!nodes %in% (nodes-1)]
# if(2 %in% nodes){
# nodes[which(nodes == 2)] = 1
# }
# important.nodes = nodes
#
# n.xylon = xy_smoothed[[1]][,1]
# n.xylat = xy_smoothed[[1]][,2]
# if(!1%in%important.nodes){
# n.xylon = c(N$utm.x[1], n.xylon)
# n.xylat = c(N$utm.y[1], n.xylat)
# important.nodes = c(1, important.nodes)
# }
# if(!nrow(N)%in%important.nodes){
# n.xylon = c(n.xylon, N$utm.x[nrow(N)])
# n.xylat = c(n.xylat, N$utm.y[nrow(N)])
# important.nodes = c(important.nodes, nrow(N))
# }
#
# for(i in 1:(length(important.nodes)-1)){
#
# # from = which(paste(round(N$lon, 10), round(N$lat, 10)) == paste(round(n.xylon[i],10), round(n.xylat[i], 10)))
# # to = which(paste(round(N$lon, 10), round(N$lat, 10)) == paste(round(n.xylon[i+1],10), round(n.xylat[i+1], 10)))
# from = important.nodes[i]
# to = important.nodes[i+1]
# fill.n = to-from
#
# tes = cbind(c(N$utm.y[important.nodes[i]],N$utm.y[important.nodes[i+1]]) , c(N$utm.x[important.nodes[i]], N$utm.x[important.nodes[i+1]]))
# tes.seg = st_cast( st_multipoint(tes, dim = "XY"), 'LINESTRING' )
# to.fill = st_sample(tes.seg, size = fill.n, type = "regular")
# N$new.lon[(from+1):to] = to.fill[[1]][,1]
# N$new.lat[(from+1):to] = to.fill[[1]][,2]
# }
# N$new.lat[important.nodes] = n.xylon
# N$new.lon[important.nodes] = n.xylat
#
#
# xyt = st_cast( st_multipoint(cbind(N$new.lat, N$new.lon), dim = "XY"), 'LINESTRING' )
# cxyllt = st_sfc(xyt, crs = "+proj=utm +zone=20 +ellps=GRS80 +datum=NAD83 +units=m" )
# herex = st_transform(cxyllt, CRS("+init=epsg:4326"))
#
#
# N$old.lat = N$lat
# N$old.lon = N$lon
# N$lat = herex[[1]][,2]
# N$lon = herex[[1]][,1]
N$lat = N$nlat
N$lon = N$nlong
}
if(nrow(N)>1) {
# t1 gives the approximate time of net lift-off from bottom
# this is not good enough as there is a potential backdrift period before net lift off
# this means distance trawled calculations must use the geo-positioning of the boat
# at maximum tension and not the position at net movemnent up
# (otherwise, this would introduce a net bias towards smaller swept areas).
pos = c( "lon", "lat" )
distance.from.start = geosphere::distCosine(N[1, pos], N[, pos])/1000 # in km^2
end = which.max( distance.from.start)
if( !is.finite(end) ) end=nrow(N)
n = N[ 1:end , ]
out$vel = mean(n$speed, na.rm=T, trim=0.1)
out$vel_sd = sd(n$speed, na.rm=T)
out$slon=n$lon[1]
out$slat=n$lat[1]
out$netmind_n=end
delta.distance = NULL
n$distance = NA
if (nrow(n) > 10) {
# integrate area:: piece-wise integration is used as there is curvature of the fishing track (just in case)
delta.distance = geosphere::distGeo ( n[ 1:(end-1), pos ], n[ 2:end, pos ] ) / 1000 ## in meters convert to km
n$distances = c( 0, cumsum( delta.distance ) ) # cumsum used to do piecewise integration of distance
# model/smooth/interpolate the spreads
n$doorspread.predicted = NA
ii = which( is.finite( n$doorspread ) )
if ( length(ii) > n.req ) {
# recall that doorspread is actually wingspread ...
nn = as.data.table( n[ii,] )
nn = nn[, median(doorspread, na.rm=TRUE), by=(distances)]
n$doorspread.predicted = approx( x=nn$distances, y=nn$V1, xout=n$distances, method="linear", rule=2, na.rm=TRUE )$y
#turned off gam model in December 20, 2013 giving unrealistic values for spread as the
# new esnoar files have 0 and NA whereas older netmind are filled with previous value
#gam.model = try( gam( doorspread ~ s(distances, k=5, bs="ts"), data=n[ii,], optimizer=c("outer", "nlm")), silent = T )
#if ( ! "try-error" %in% class( gam.model )) {
# n$doorspread.predicted = predict( gam.model, newdata=n, newdata.guaranteed=T )
#} else {
# n$doorspread.predicted = approx( x=n$distances, y=n$doorspread, xout=n$distances, method="linear", rule=2 )$y
#}
}
if ( length( which( is.finite( n$doorspread.predicted ) ) ) < 10 ) {
n$doorspread.predicted = mean( n$doorspread , na.rm=T, trim=0.1 )
}
mean.doorspreads = ( n$doorspread.predicted[1:(end-1)] + n$doorspread.predicted[2:(end)] ) / 2 / 1000 # mean between two ends
partial.area = delta.distance * mean.doorspreads
out$surfacearea = sum( partial.area ) # km^2
out$surfacearea = abs( out$surfacearea )
out$spread = mean(n$doorspread.predicted, na.rm=T, trim=0.1)/1000 # in km
spread_sd = sd(n$doorspread.predicted, na.rm=T )/1000
if(!is.na(spread_sd) & spread_sd!=0) out$spread_sd = spread_sd #if just using the mean from above do not over write spread_sd
out$distance=n$distances[end]
# 2022 Temperature now available on Marport Sensors. Keep last 3/4 of bottom contact to help parse out latency
temp.to.keep = n$temperature[(length(n$temperature)/4):length(n$temperature)]
if ( length(which(!is.na(temp.to.keep))) > 4){
out$temperature.n = mean(temp.to.keep, na.rm=T, trim=0.1)
out$temperature_sd.n = sd(temp.to.keep, na.rm=T)
}
else{
out$temperature.n = NA
out$temperature_sd.n = NA
}
}
}
print(out)
return (out)
}
jae0/bio.snowcrab documentation built on Nov. 6, 2024, 10:10 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 net.configuration
net.configuration = function( N, t0=NULL, t1=NULL, set_timestamp=NULL, yr=NULL, plotdata=TRUE ) {
# N is netmind data
# t0 is current best estimate of start and end time
# set_timestamp is timestamp from set-log .. alternate if there is no other useful time marker
# note:: NETMIND data stream calls the wingspread as "DOORSPREAD".
# It is not. It is wingspread.
plotdir = project.datadirectory("bio.snowcrab", "data", "netmind", "figures" )
# create default output should the following fail
out = data.frame( slon=NA, slat=NA, distance=NA, spread=NA, spread_sd=NA,
surfacearea=NA, vel=NA, vel_sd=NA, netmind_n=NA, t0=NA, t1=NA, dt=NA, yr=NA,
temperature.n=NA, temperature_sd.n=NA )
n.req = 10
if (!is.null(t0)) {
if (length(t0)>1) {t0 = NULL}
if (is.na(t0)) t0 = NULL
}
if (!is.null(t1)) {
if (is.na(t1)) t1 = NULL
}
#changed this switch from depgth filed to lat as if there is not depth info can still run script
if ( length( which( is.finite( N$lat))) < n.req ) print(N[1,])
problem = F
# time checks
if ( is.null(t0) ){
if ( !is.null(set_timestamp) ) {
t0 = set_timestamp # no data in t0 ,, use set_timestamp as alternate
set_timestamp =NULL # remove to cause no more effects
}
}
#bad.list = c('netmind.S26092014.9.541.17.48.304',
# 'netmind.S20092007.8.333.17.27.241' )
bad.list = NULL
bad.list = unique( c(bad.list, p$netmensuration.problems) )
if ( N$netmind_uid[1] %in% bad.list) {
return(out)
}
if ( any( is.null( t1 ) || is.null(t0) ) ) {
# try to determine from netmind data if minilog/seadbird data methods have failed. .. not effective due to noise/and small data stream
print(N$netmind_uid[1])
M = N[, c("timestamp", "depth") ]
oo = which( is.finite( M$depth ))
if ( length(oo) < n.req ) return(out)
if(!is.null(set_timestamp)) settimestamp=set_timestamp
if(is.null(set_timestamp)) settimestamp=t0
time.gate = list( t0=settimestamp - dminutes(6), t1=settimestamp + dminutes(12) )
bcp = list(id=N$netmind_uid[1], nr=nrow(M), YR=yr, trip = unlist(strsplit(N$netmind_uid[1], "\\."))[2], tdif.min=3, tdif.max=11, time.gate=time.gate,
depth.min=20, depth.range=c(-25,15), eps.depth=3,
smooth.windowsize=5, modal.windowsize=5, datasource = "snowcrab",
noisefilter.trim=0.05, noisefilter.target.r2=0.8, noisefilter.quants=c(0.05, 0.95) )
if(yr<2007)bcp$from.manual.archive=FALSE # manual touchdown only done since 2007
bcp = bottom.contact.parameters( bcp ) # add other default parameters
#BC: Determine if this station was done yet, if not then we want user interaction.
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) bcp$user.interaction = FALSE
else bcp$user.interaction = TRUE
}
dp.out.range = which(M$depth < .1 | M$depth > 360)
if(length(dp.out.range)>0){
M$depth[dp.out.range] = NA
}
M$wingspread = M$doorspread
bc = NULL
bc = bottom.contact( x=M, bcp=bcp )
if ( is.null(bc) || (!is.null( bc$res) && ( ( !is.finite(bc$res$t0 ) || !is.finite(bc$res$t1 ) ) ) )) {
bcp$noisefilter.target.r2=0.8
bc = bottom.contact( x=M, bcp=bcp )
}
if ( is.null(bc) || (!is.null( bc$res) && ( ( !is.finite(bc$res$t0 ) || !is.finite(bc$res$t1 ) ) ) )) {
# try once more with random settings
bcp$noisefilter.inla.h = 0.05
bcp$noisefilter.target.r2= 0.75
bc = bottom.contact( x=M, bcp=bcp )
}
if ( is.null(bc) || (!is.null( bc$res) && ( ( !is.finite(bc$res$t0 ) || !is.finite(bc$res$t1 ) ) ) )) {
# try once more with more random settings .. noise is high in netminds data
bcp$eps.depth = 5
M$depth = jitter( M$depth, amount = bcp$eps.depth/10 )
bcp$noisefilter.inla.h = 0.1
bc = bottom.contact( x=M, bcp=bcp )
}
if ( !is.null(bc) ) {
if (plotdata) {
bottom.contact.plot( bc )
plotfn = file.path( plotdir, paste(bcp$id, "pdf", sep="." ) )
print (plotfn)
dev.flush()
dev.copy2pdf( file=plotfn )
graphics.off()
}
}
if (is.null(t0) & !is.null(bc$bottom0) ) t0 = bc$bottom0
if (is.null(t1) & !is.null(bc$bottom1) ) t1 = bc$bottom1
N = N[ bc$bottom.contact , ]
}
if (all(is.na(t0))) t0=NA
if (all(is.na(t1))) t1=NA
if ( is.null(t1) || !is.finite(t1) ) {
t1_tmp = NA # rather than guessing, flag and then fill later
} else {
t1_tmp = t1
}
# if we are here, it is because a reasonable amount of data is present ..
# do some more checks and get a first estimate of some parameters in case other errors/limits are found
out$slon=N$lon[1]
out$slat=N$lat[1]
out$spread=mean( N$doorspread, na.rm=T ) / 1000 # though called doorspread it is actually wingspread
out$spread_sd=sd( N$doorspread, na.rm=T ) /1000
# first set the impossible door spreads to NA and then interpolate based upon a 5 and 95 % quantiles
# Netmind has very noisy data
hl = c( 3, 16 ) # m
ihl = which( N$doorspread < hl[1] | N$doorspread > hl[2] )
if (length (ihl)>0 ) N$doorspread[ihl] = NA
quantiles.to.trim = c(0.05, 0.95)
qnt = quantile( N$doorspread[N$doorspread>0], quantiles.to.trim, na.rm=T)
iqnt = which( N$doorspread < qnt[1] | N$doorspread > qnt[2] )
N$doorspread[ iqnt ] = NA
gooddoor = which( is.finite(N$doorspread) )
baddoor = which( ! is.finite(N$doorspread) )
if (length(baddoor)>0) N$doorspread[ baddoor ] = NA
if ( length( gooddoor) < n.req ) {
#problem =T turned this off December 23, 2013 as if we have decent gps data now we can calc distance
if ( length( gooddoor) >0 ) {
out$netmind_n <- length(gooddoor)
out$spread = mean( N$doorspread[ gooddoor ], na.rm=T ) / 1000
out$spread_sd = sd( N$doorspread[ gooddoor ], na.rm=T ) / 1000
}
}
if ( length(t0)>1) t0 = NA
out$t0 = t0
out$t1 = t1
out$dt = difftime( as.POSIXct(t1), as.POSIXct(t0) ) # minutes
if(!is.na(out$t0)) out$yr = lubridate::year( out$t0)
if(is.na(out$t0)) out$yr = lubridate::year(N$timestamp[1])
itime = which( N$timestamp >= t0 & N$timestamp <= t1 )
if ( length( itime) < n.req ) problem = T
if (problem) {
## Need to check on these at some point
out$t0 = t0
out$t1 = NA
out$dt = NA
return(out)
}
# eOR checks
# Higher GPS data accuracy in 2019 resulted in high approximation of distance travelled. This was due to frequent occasions of backwards
# motion (assuming from boat roll and pitch) which was added to the total distance. When looking at previous years this also seem to have
# but not to the same extent however removing the backwards motion would result in better data accuracy from year using e-sonar 2014+.
# The solution is to trim the data by 3/4, or, with the current recording rate of every 4 seconds instead of every 1 second.
# The above works when the metrics are repeated but not when we actually only record the values when we get a sensor hit. This does not
# work for the 2021 data. Trying to smooth latitude and longitude so that we dont have data loss.
# This again did not work as paths were way out in some cases. Resolved by applying a windowed average smoothing
# Keeping code for st_simplify in case
itime.smoothing = which( N$timestamp >= (t0 - lubridate::seconds(10)) & N$timestamp <= (t1 + lubridate::seconds(10) ))
rolling.ave = 6
Ntemp = N[ itime.smoothing,]
Ntemp$nlat = NA
Ntemp$nlong = NA
for(i in 1:nrow(Ntemp)){
inds = (i-rolling.ave):(i+rolling.ave)
inds = inds[which(inds>0)]
inds = inds[which(inds<=nrow(Ntemp))]
Ntemp$nlong[i] = mean(Ntemp$lon[inds])
Ntemp$nlat[i] = mean(Ntemp$lat[inds])
}
#plot(Ntemp$lon,Ntemp$lat, type = "l", col = "#000000E6", lwd = 3)
#lines(Ntemp$nlon,Ntemp$nlat, col = 'red', lwd = 2)
itime = which( Ntemp$timestamp >= t0 & Ntemp$timestamp <= t1 )
N = Ntemp[ itime, ]
if(yr>=2014){
# require(sf)
# tolerence = 5
# N$new.lat = NA
# N$new.lon = NA
# xy = st_cast( st_multipoint(cbind(N$lon,N$lat), dim = "XY"), 'LINESTRING' )
# cxyll = st_sfc(xy, crs = "+init=epsg:4326" )
# c_xy = st_transform(cxyll, st_crs("+proj=utm +zone=20 +ellps=GRS80 +datum=NAD83 +units=m" ))
# N$utm.y = c_xy[[1]][,2]
# N$utm.x = c_xy[[1]][,1]
# xy_smoothed = st_simplify( c_xy, dTolerance=tolerence)
# tol.count = tolerence+1
# while(any(grepl("MULTILINESTRING", class(xy_smoothed)))){
# tol.count = tol.count+1
# print(paste("Moving tolerence to:", tol.count))
# xy_smoothed = st_simplify( c_xy, dTolerance=tol.count)
# if(tol.count > 100) break()
# }
#
#
#
# nodes = which(paste(N$utm.y, N$utm.x) %in% paste(xy_smoothed[[1]][,2],xy_smoothed[[1]][,1]))
# nodes = nodes[!nodes %in% (nodes-1)]
# if(2 %in% nodes){
# nodes[which(nodes == 2)] = 1
# }
# important.nodes = nodes
#
# n.xylon = xy_smoothed[[1]][,1]
# n.xylat = xy_smoothed[[1]][,2]
# if(!1%in%important.nodes){
# n.xylon = c(N$utm.x[1], n.xylon)
# n.xylat = c(N$utm.y[1], n.xylat)
# important.nodes = c(1, important.nodes)
# }
# if(!nrow(N)%in%important.nodes){
# n.xylon = c(n.xylon, N$utm.x[nrow(N)])
# n.xylat = c(n.xylat, N$utm.y[nrow(N)])
# important.nodes = c(important.nodes, nrow(N))
# }
#
# for(i in 1:(length(important.nodes)-1)){
#
# # from = which(paste(round(N$lon, 10), round(N$lat, 10)) == paste(round(n.xylon[i],10), round(n.xylat[i], 10)))
# # to = which(paste(round(N$lon, 10), round(N$lat, 10)) == paste(round(n.xylon[i+1],10), round(n.xylat[i+1], 10)))
# from = important.nodes[i]
# to = important.nodes[i+1]
# fill.n = to-from
#
# tes = cbind(c(N$utm.y[important.nodes[i]],N$utm.y[important.nodes[i+1]]) , c(N$utm.x[important.nodes[i]], N$utm.x[important.nodes[i+1]]))
# tes.seg = st_cast( st_multipoint(tes, dim = "XY"), 'LINESTRING' )
# to.fill = st_sample(tes.seg, size = fill.n, type = "regular")
# N$new.lon[(from+1):to] = to.fill[[1]][,1]
# N$new.lat[(from+1):to] = to.fill[[1]][,2]
# }
# N$new.lat[important.nodes] = n.xylon
# N$new.lon[important.nodes] = n.xylat
#
#
# xyt = st_cast( st_multipoint(cbind(N$new.lat, N$new.lon), dim = "XY"), 'LINESTRING' )
# cxyllt = st_sfc(xyt, crs = "+proj=utm +zone=20 +ellps=GRS80 +datum=NAD83 +units=m" )
# herex = st_transform(cxyllt, CRS("+init=epsg:4326"))
#
#
# N$old.lat = N$lat
# N$old.lon = N$lon
# N$lat = herex[[1]][,2]
# N$lon = herex[[1]][,1]
N$lat = N$nlat
N$lon = N$nlong
}
if(nrow(N)>1) {
# t1 gives the approximate time of net lift-off from bottom
# this is not good enough as there is a potential backdrift period before net lift off
# this means distance trawled calculations must use the geo-positioning of the boat
# at maximum tension and not the position at net movemnent up
# (otherwise, this would introduce a net bias towards smaller swept areas).
pos = c( "lon", "lat" )
distance.from.start = geosphere::distCosine(N[1, pos], N[, pos])/1000 # in km^2
end = which.max( distance.from.start)
if( !is.finite(end) ) end=nrow(N)
n = N[ 1:end , ]
out$vel = mean(n$speed, na.rm=T, trim=0.1)
out$vel_sd = sd(n$speed, na.rm=T)
out$slon=n$lon[1]
out$slat=n$lat[1]
out$netmind_n=end
delta.distance = NULL
n$distance = NA
if (nrow(n) > 10) {
# integrate area:: piece-wise integration is used as there is curvature of the fishing track (just in case)
delta.distance = geosphere::distGeo ( n[ 1:(end-1), pos ], n[ 2:end, pos ] ) / 1000 ## in meters convert to km
n$distances = c( 0, cumsum( delta.distance ) ) # cumsum used to do piecewise integration of distance
# model/smooth/interpolate the spreads
n$doorspread.predicted = NA
ii = which( is.finite( n$doorspread ) )
if ( length(ii) > n.req ) {
# recall that doorspread is actually wingspread ...
nn = as.data.table( n[ii,] )
nn = nn[, median(doorspread, na.rm=TRUE), by=(distances)]
n$doorspread.predicted = approx( x=nn$distances, y=nn$V1, xout=n$distances, method="linear", rule=2, na.rm=TRUE )$y
#turned off gam model in December 20, 2013 giving unrealistic values for spread as the
# new esnoar files have 0 and NA whereas older netmind are filled with previous value
#gam.model = try( gam( doorspread ~ s(distances, k=5, bs="ts"), data=n[ii,], optimizer=c("outer", "nlm")), silent = T )
#if ( ! "try-error" %in% class( gam.model )) {
# n$doorspread.predicted = predict( gam.model, newdata=n, newdata.guaranteed=T )
#} else {
# n$doorspread.predicted = approx( x=n$distances, y=n$doorspread, xout=n$distances, method="linear", rule=2 )$y
#}
}
if ( length( which( is.finite( n$doorspread.predicted ) ) ) < 10 ) {
n$doorspread.predicted = mean( n$doorspread , na.rm=T, trim=0.1 )
}
mean.doorspreads = ( n$doorspread.predicted[1:(end-1)] + n$doorspread.predicted[2:(end)] ) / 2 / 1000 # mean between two ends
partial.area = delta.distance * mean.doorspreads
out$surfacearea = sum( partial.area ) # km^2
out$surfacearea = abs( out$surfacearea )
out$spread = mean(n$doorspread.predicted, na.rm=T, trim=0.1)/1000 # in km
spread_sd = sd(n$doorspread.predicted, na.rm=T )/1000
if(!is.na(spread_sd) & spread_sd!=0) out$spread_sd = spread_sd #if just using the mean from above do not over write spread_sd
out$distance=n$distances[end]
# 2022 Temperature now available on Marport Sensors. Keep last 3/4 of bottom contact to help parse out latency
temp.to.keep = n$temperature[(length(n$temperature)/4):length(n$temperature)]
if ( length(which(!is.na(temp.to.keep))) > 4){
out$temperature.n = mean(temp.to.keep, na.rm=T, trim=0.1)
out$temperature_sd.n = sd(temp.to.keep, na.rm=T)
}
else{
out$temperature.n = NA
out$temperature_sd.n = NA
}
}
}
print(out)
return (out)
}
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.