R/calculateBiomass_Simmonds2010.R

In imarpe/path: Pelagic Assessment Tools for Hydroacoustics

# Calculate biomass using Simmonds et al. 2010 method ---------------------
# TO_DO: Brief explanation
# Main method, auxiliar functions after this definition
.calculateBiomass_simmonds2010 = function(data, sp, boot, n, verbose,
                                          coverage, dateFormat, ...) {
  
  if(!isTRUE(boot)) n = 1
  
  specie = getSpeciesInfo(sp) # get species information
  marcas = .createMarks(specie)
  
  # indCols = c("interval", "lon", "lat", "area", "fecha", specie$acustica)
  baseAcustica = data$bitacoraAcustica #[, indCols]
  baseCalas    = data$baseCalas
  
  coverage = .coverageIndex(baseAcustica, sp, method=coverage)

  # Calcular los valores de frecuencia ponderada
  frecuenciaPonderada = .getLength(baseCalas = baseCalas, sp = sp)
  TS = .calculateTS(sp)
  
  infoAreasThis = merge(data.frame(code = unique(baseAcustica$area)), 
                        infoAreasIsop, sort = FALSE)
  listaAreas = infoAreasThis$code
  areaMN     = infoAreasThis$area_mn2
  
  # Definir bases de datos finales
  
  ecoResampleBase = matrix(nrow = length(listaAreas), ncol = n,
                           dimnames=list(listaAreas, seq_len(n)))
  
  biomassResampleBase = matrix(nrow = length(listaAreas), ncol = n,
                               dimnames=list(listaAreas, seq_len(n)))
  
  abundanceResampleBase = matrix(nrow = length(listaAreas), ncol = n,
                                 dimnames=list(listaAreas, seq_len(n)))
  
  bootMatrixAreaN = array(NA, dim = c(length(listaAreas), length(marcas), n))
  bootMatrixAreaB = array(NA, dim = c(length(listaAreas), length(marcas), n))
  
  datesBYarea = data.frame(date = as.Date(data$bitacoraAcustica$date_m, format = dateFormat),
                           area = data$bitacoraAcustica$area)
  datesBYarea = as.character(as.Date(with(datesBYarea, tapply(date, area, mean)), 
                                     origin = "1970-1-1"))
  
  
  if(isTRUE(boot)) cat("\nStarting bootstraping. \t n = ", n, "\n\n")
  
  for(j in seq_len(n)) {
    
    if(isTRUE(verbose) & isTRUE(boot)) 
      cat("Iteration", j, "of", n, "\n")
    
    baseTallasN = mat.or.vec(nr = length(listaAreas), nc = length(marcas))
    baseTallasB = mat.or.vec(nr = length(listaAreas), nc = length(marcas))
    baseAreaBiomasa     = numeric(length(listaAreas))
    baseAreaAbundancia  = numeric(length(listaAreas))
    ecoabundance        = numeric(length(listaAreas))
    
    # Estimacion de B, N, frecuencia de tallas y otros por area
    for (i in seq_along(listaAreas)) {
      
      # area i information
      thisArea               = listaAreas[i] 
      baseAcustica.thisArea  = baseAcustica[baseAcustica$area == thisArea,]
      baseAcustica.thisArea  = baseAcustica.thisArea[!is.na(baseAcustica.thisArea$areaCercana),]
      calas.thisArea         = which(baseCalas$area == thisArea)
      coverage.thisArea      = coverage[as.character(thisArea)]
      
      # hauls for this area
      lances = .getLances(baseAcustica.thisArea, baseCalas, calas.thisArea, sp)
      baseCalas.thisArea     = baseCalas[lances, ]
      
      # resampling
      newLances       = .resampleLances(lances, boot=boot)
      newBaseCalas    = baseCalas[newLances, ]
      newBaseAcustica = .resampleTracks(baseAcustica.thisArea, boot=boot)
      newBaseAcustica = .resampleProportions(newBaseAcustica, newBaseCalas, 
                                             sp=specie$acustica, boot=boot)
      
      # mean length frequency and ecoabundance
      vectorFrecuencias = .calculateRelativeAbundance(frecuenciaPonderada, newLances)
      meanEcoabundance  = mean(newBaseAcustica[, specie$acustica]) # mean ecoabundance
      
      # Abundance estimation (by length)
      sigma         = 4*pi*sum(vectorFrecuencias*(10^(TS/10)))
      abundancia    = (meanEcoabundance/sigma)*areaMN[i]*coverage.thisArea
      tallasNumero  = vectorFrecuencias*abundancia*1e-6
      weights       = .getWeights(baseCalas, lances, sp)            
      tallasBiomasa = weights*tallasNumero  
      biomasa       = sum(tallasBiomasa, na.rm = TRUE) # Biomasa total
      
      # save results
      baseTallasN[i,]       = tallasNumero
      baseTallasB[i,]       = tallasBiomasa
      baseAreaBiomasa[i]    = biomasa
      baseAreaAbundancia[i] = sum(abundancia, na.rm = TRUE)
      ecoabundance[i]       = meanEcoabundance
      
    } # end of listaAreas loop
    
    ecoResampleBase[, j]       = ecoabundance
    biomassResampleBase[, j]   = baseAreaBiomasa
    abundanceResampleBase[, j] = baseAreaAbundancia
    bootMatrixAreaN[, , j]     = baseTallasN
    bootMatrixAreaB[, , j]     = baseTallasB
    
  } # end of bootstrap
  
  # Preparacion de las bases de datos
  
  bootData   = NULL
  
  if(isTRUE(boot)) {
    
    # Armar salidas
    bootData = list(biomassResampleBase = biomassResampleBase,
                    abundanceResampleBase = abundanceResampleBase,
                    ecoResampleBase = ecoResampleBase,
                    bootMatrixAreaN = bootMatrixAreaN, 
                    bootMatrixAreaB = bootMatrixAreaB,
                    biomass = colSums(biomassResampleBase),
                    listaAreas = listaAreas)
    
    baseTallasN        = apply(bootMatrixAreaN, c(1,2), mean, na.rm=TRUE)
    baseTallasB        = apply(bootMatrixAreaB, c(1,2), mean, na.rm=TRUE)
    baseAreaBiomasa    = rowMeans(biomassResampleBase, na.rm=TRUE)
    baseAreaAbundancia = rowMeans(abundanceResampleBase, na.rm=TRUE)
    ecoabundance       = rowMeans(ecoResampleBase, na.rm=TRUE)
    
  } 
  
  # remove NaN, check for a neat solution
  baseTallasN[is.na(baseTallasN)] = 0
  baseTallasB[is.na(baseTallasB)] = 0
  
  # base por areas
  baseArea = merge(data.frame(code=listaAreas), infoAreasThis, sort=FALSE)
  baseArea$nasc       = ecoabundance 
  baseArea$biomass    = baseAreaBiomasa
  baseArea$abundance  = baseAreaAbundancia
  baseArea$densityN   = baseArea$abundance/baseArea$area_mn2
  baseArea$densityB   = baseArea$biomass/baseArea$area_mn2
  baseArea$date       = datesBYarea[match(baseArea$code, names(datesBYarea))]
  
  # base por Areas y Tallas
  output = merge(data.frame(code=rep(listaAreas, each=length(marcas))), 
                 infoAreasThis, sort=FALSE)
  output = data.frame(output, 
                      L=marcas, 
                      abundance=as.numeric(t(baseTallasN)), 
                      biomass=as.numeric(t(baseTallasB)))    
  output$densityN = output$abundance/output$area_mn2
  output$densityB = output$biomass/output$area_mn2
  output$date     = datesBYarea[match(output$code, names(datesBYarea))]
  output$coverage = coverage
  
  # TO_DO (general output path format ^.^): 
  # only raw as outputs, final output in estimarBiomasa
  # Asignar nombres de columna
  colnames(baseTallasN) = marcas
  colnames(baseTallasB) = marcas
  
  raw = list(base=output, area=baseArea, 
             baseTallasN = baseTallasN, baseTallasB = baseTallasB, 
             boot=bootData)
  
  out = list(raw = raw, results = .groupData(output))
  
  return(out)
}


# Auxiliar function for getBiomass_Simmonds2010 ---------------------------

.getTracks = function(base, dtLim = 10, dxLim = 2) {
  # Definir pasos de tiempo y espacio entre un punto y otro
  # TO_DO: nrow(base)<2, use diff.
  dt = base$fecha[2:nrow(base)] - base$fecha[1:(nrow(base) - 1)]
  dx = .distancia(x1 = base$lon[2:nrow(base)], 
                  x2 = base$lon[1:(nrow(base) - 1)], 
                  y1 = base$lat[2:nrow(base)], 
                  y2 = base$lat[1:(nrow(base) - 1)])
  
  # Definir los giros: Puntos en los que dt > dtLim & dx > dxLim
  turns = which(dt > dtLim & dx > dxLim) # antes estaba como dx > dtLim
  turns = c(0, turns, nrow(base))
  tracks = rep(seq_along(turns[-1]), diff(turns))
  return(tracks)
}

.resampleTracks = function(base, boot) {
  
  if(!isTRUE(boot)) return(base)
  
  tracks  = .getTracks(base)
  
  if(length(unique(tracks))==1) return(base)
  
  muestra = sample(x = as.numeric(unique(tracks)), size = length(unique(tracks)), 
                   replace = TRUE, prob = as.numeric(table(tracks)/length(tracks)))
  
  newBase = NULL
  
  for(k in seq_along(muestra)) {
    newBase = rbind(newBase, base[tracks == muestra[k], ])
  }
  
  return(newBase)
  
}

.resampleProportions = function(base, calas, sp, boot) {
  
  if(!isTRUE(boot)) return(base)  
  
  nasc      = base[, sp]
  totalNasc = nasc/base$prop
  newNasc   = totalNasc*mean(calas$prop, na.rm=TRUE)
  newNasc[!base$isAllocated] = nasc[!base$isAllocated]*runif(1)/0.5
  
  base[, sp] = newNasc
  
  return(base)
}

.resampleCalas = function(calas, boot) {
  
  if(!isTRUE(boot)) return(calas)
  
  ind = sample(nrow(calas), replace=TRUE)
  
  return(calas[ind, ])
}

.resampleLances = function(lances, boot) {
  
  if(!isTRUE(boot)) return(lances)
  
  newLances = sample(lances, replace=TRUE)
  
  return(newLances)
}

.getGravityCenter = function(base, sp) {
  if(sum(base[, sp], na.rm = TRUE)==0) {
    lonCentroGravedad = mean(base$lon)
    latCentroGravedad = mean(base$lat)
  } else {
    lonCentroGravedad = sum(base$lon*base[, sp])/sum(base[, sp])
    latCentroGravedad = sum(base$lat*base[, sp])/sum(base[, sp])   
  }
  
  return(list(lat=latCentroGravedad, lon=lonCentroGravedad))
}

.calaCercana = function(acustica, calas, sp) {
  GC = .getGravityCenter(base=acustica, sp=sp)
  distancias = .distancia(x1 = GC$lon, x2 = calas$lon, 
                          y1 = GC$lat, y2 = calas$lat)
  # Identificar la posicion del punto con la menor distancia
  cercana = which.min(distancias)
  return(cercana)
}

.calculateRelativeAbundance = function(frecuenciaPonderada, lances) {
  vectorFrecuencias = if(length(lances) == 1)
    frecuenciaPonderada[lances,] else 
      colSums(frecuenciaPonderada[lances,], na.rm = TRUE)
  vectorFrecuencias = .relative(vectorFrecuencias)
  return(vectorFrecuencias)
}

.getLances = function(acustica, calas, lances, sp) {
  sp = getSpeciesInfo(sp)$acustica
  output = if(length(lances)>0) 
    lances else .calaCercana(acustica, calas, sp)
  return(output)
} 

.getWeights = function(baseCalas, lances, sp) {
  # if(sp!="anchoveta") stop("Currently working with anchovy only")
  marcas = .createMarks(getSpeciesInfo(sp))
  a = mean(baseCalas$a[lances], na.rm=TRUE)
  b = mean(baseCalas$b[lances], na.rm=TRUE)
  weights = a*marcas^b
  return(weights)
}