R/crucero_funciones.R

In imarpe/path: Pelagic Assessment Tools for Hydroacoustics

getSpeciesInfo = function(sp, data=NULL) {
  
#   otros = list(Lmin=min(data$talla-2, na.rm=TRUE), 
#                Lmax=max(data$talla+2, na.rm=TRUE), 
#                bin=1, diplay.name=sp, juvenile=NULL, unit="cm", 
#                lengthType="total")

  otros = NULL
  
  out = if(sp %in% rownames(species)) as.list(species[sp, ]) else otros
  
  return(out)
  
}

.getNewAxis = function(pos, marks) {
  newAxis = lm(pos ~ marks)
  out = function(x) {
    predict(newAxis, newdata=list(marks=x))
  }
  return(out)
}

.createMarksClasses = function(specie) {
  marks = seq(from=specie$Lmin, to=specie$Lmax, by=specie$bin)
  freqs = numeric(length(marks))
  names(freqs) = marks
  return(freqs)
}

.createMarks = function(specie, phi=FALSE) {
  marks = seq(from=specie$Lmin, to=specie$Lmax, by=specie$bin)
  if(isTRUE(phi)) {
    marks_inf = marks - 0.5*specie$bin
    marks_sup = marks + 0.5*specie$bin
    marks = sort(unique(c(marks_inf, marks_sup)))
  }
  return(marks)
}

.removeTails = function(x) {
  y = sign(abs(x))
  r = rle(y)
  if(r$value[1]==0 & r$length[1]>1) x = tail(x, -r$length[1]+1)
  if(tail(r$value, 1)==0 & tail(r$length, 1)>1) x = head(x, -tail(r$length, 1)+1)
  return(x)
}

.getModes = function(x, thr=0.01) {

  x = .removeTails(x)
  y = rle(x)
  yValue  = y$value
  yLength = which(y$length!=1)
  values = as.numeric(names(yValue))
  newValues = rowMeans(cbind(values[yLength-1], values[yLength+1]))
  values[yLength] = newValues
  names(yValue) = values
  xdif  = rle(sign(diff(yValue)))$value
  modes = as.numeric(names(xdif[xdif==1]))
  x = x/sum(x, na.rm=TRUE)
  modes = c(na.omit(modes[x[as.character(modes)]>thr]))
  return(modes)
  
}

plotFreq = function(sp, data, ...) {
  
  x = data[data$especie==sp, ] # filter data
  
  specie = getSpeciesInfo(sp, data=x) # get species info
  y = tapply(x$freq, INDEX = x$talla, FUN=sum) # total freqs!
  imarks = as.numeric(names(y)) # observed marks of class
  
  freqs = .createMarksClasses(specie) # create full classes
  marks = .createMarks(specie)
  freqs[marks %in% imarks] = y # assign observed values
  n = sum(freqs, na.rm=TRUE)
  ran = range(imarks)
  juv = 100*sum(freqs[marks<=specie$juvenile], na.rm=TRUE)/n
  freqs = 100*freqs/sum(freqs, na.rm=TRUE)
  modes = .getModes(freqs)
  modesText = paste(modes, collapse=", ")
  ylim = c(0, 1.25*max(freqs, na.rm=TRUE))
  pos = barplot(freqs, ylim=ylim, axes=FALSE, axisname=FALSE,
                border=NA, space=0.1, ...)
  newAxis = .getNewAxis(pos, marks)
  axis(1, at=pos, labels = marks)
  axis(2, las=2)
  # axis labels
  mtext("Porcentaje (%)", 2, line=2.5, cex=1.2)
  mtext(sprintf("Longitud (%s, %s)", specie$lengthType, specie$unit),
        1, line=2.5, cex=1.2)
  # species info
  mtext(specie$diplay.name, 3, line=-2.0, adj=0.05, cex=1.2)
  mtext(specie$sci.name, 3, line=-3.2, adj=0.05, cex=1.0, font=3)
  # data info
  mtext(sprintf("n = %d", n), 3, line=-2.0, adj=0.95)
  mtext(sprintf("rango = %d - %d %s", ran[1], ran[2], specie$unit), 
        3, line=-3.2, adj=0.95)
  if(!is.null(specie$juvenile)) 
    mtext(sprintf("%.01f%% juveniles", juv), 3, line=-4.4, adj=0.95)
  mtext("Modas:", 3, line=-5.7, adj=0.96)
  mtext(sprintf("%s %s", modesText, specie$unit), 3, line=-7.0, adj=0.95)
  if(!is.null(specie$juvenile)) 
    abline(v=newAxis(specie$juvenile), lty=2, col="red")
  box()
  
  return(invisible())
}

calculateConditionFactor = function(data, sp, ...) {

  if(is.null(data$baseBiologico) & is.null(data$baseBiometrico))
    stop("Cannot calculate condition factor from this dataset, use average values")
  
  ab_biol = .getConditionFactor(data=data$baseBiologico, sp=sp)
  ab_biom = .getConditionFactor(data=data$baseBiometrico, sp=sp)
  ab_base = .getConditionFactor(sp=sp)  # TO_DO: check best value for ab_base
  
  if(!is.null(data$baseBiologico) & !is.null(data$baseBiometrico)) {
    ab = merge(ab_biol, ab_biom, sort=FALSE)
    a = ab$a.x
    a[is.na(a)] = ab$a.y[is.na(a)]
    a[is.na(a)] = ab_base$a
    
    b = ab$b.x
    b[is.na(b)] = ab$b.y[is.na(b)]
    b[is.na(b)] = ab_base$b
    
    output = data.frame(buque=ab$buque, lance=ab$lance, a=a, b=b)
    return(output)
    
  } 

  ab = if(!is.null(data$baseBiologico)) ab_biol else ab_biom

    ab$a[is.na(ab$a)] = ab_base$a
    ab$b[is.na(ab$b)] = ab_base$b

  return(ab)

}


.getConditionFactor = function(data = NULL, sp, ...) {
  
  dataName = deparse(substitute(data))
  
  if(is.null(data)) {
    ab = getSpeciesInfo(sp)[c("a", "b")]
    return(ab)
  }

  check = c("weight", "lance", "length", "buque") %in% names(data)
  if(!all(check)) stop(sprintf("Check variables in base %s", dataName))
  
  if(all(is.na(data$weight))) {
    ab = getSpeciesInfo(sp)[c("a", "b")]
    ab = data.frame(buque=data$buque, lance=data$lance, a=ab$a, b = ab$b)
    return(ab)
  }
  
  data = data[data$sp == sp, ]
  data$buqueLance = paste(data$buque, data$lance, sep="-")
  data$buqueLance = as.factor(data$buqueLance)
  model = lm(log(weight) ~ buqueLance + log(length) + 0, data=data)
  
  pars = summary(model)$coefficients[, 1:2]
  b = pars["log(length)", ]
  a = pars[grep(x=rownames(pars), pat="buqueLance"),]
  buqueLance =   do.call(rbind, strsplit(gsub(x=rownames(a), 
                                               pat="buqueLance", rep=""), 
                           split = "-"))
  ab = data.frame(buque=buqueLance[,1], lance=buqueLance[,2], 
                  a=exp(a[,1] + 0.5*a[,2]), b = b[1])
  
  return(ab)
}