R/imarpeSurvey-auxiliar.R

In imarpe/path: Pelagic Assessment Tools for Hydroacoustics

# add allocation information to bitacoraAcustica
.isAllocated = function(data) {
  bitacoraAcustica = data$bitacoraAcustica
  xdist = kali:::getDistance(data = bitacoraAcustica, ref=data$baseCalas,
                      index.return=TRUE)
  bitacoraAcustica$dist2haul = xdist[, "dist"]/1.852
  bitacoraAcustica$areaCercana = data$baseCalas$area[xdist[, "ix"]]
  bitacoraAcustica$isAllocated = bitacoraAcustica$dist2haul < 30
  return(bitacoraAcustica)
}

.addPropBySpecie = function(data, sp) {
  
  sp = getSpeciesInfo(sp)$acustica # get species information
  
  bitacoraAcustica = data$bitacoraAcustica
  propByArea = tapply(data$baseCalas$prop, INDEX = data$baseCalas$area, FUN=mean, na.rm=TRUE)
  propByArea = data.frame(areaCercana=names(propByArea), prop=propByArea)
  
  bitacoraAcustica = merge(bitacoraAcustica, propByArea, all=TRUE)
  bitacoraAcustica$spTotal = bitacoraAcustica[, sp]/bitacoraAcustica$prop
  
  return(bitacoraAcustica)
  
}


# Funcion que toma frecuencias de tallas y coeficientes para devolver un vector de 
# tallas en peso
.ponderaTallas = function(frecuencia, marcas, a, b, peso) {

  # frequence to one unit of weigth (given by 'a' units)
  output = frecuencia/sum(a*(marcas^b)*frecuencia, na.rm = TRUE)
  output = peso*output # why???????
  return(as.numeric(output))
  
}

# Criterio 1 de proporcion de cobertura: Se asumen todos con cobertura total (valor 1)
.coverageIndex_1 = function(bitacoraAcustica, sp) {
  sp = getSpeciesInfo(sp)$acustica
  data = bitacoraAcustica[, c("area", sp)]

  data = data[complete.cases(data),]
  
  acusticaPositiva  = tapply(X=data[, sp], INDEX=data$area, FUN=sum, na.rm=TRUE)
  acusticaPositiva  = names(acusticaPositiva[acusticaPositiva > 0])
  
  out = rep(1, len=length(unique(data$area)))
  names(out) = unique(data$area)
#   out[unique(data$area) %in% acusticaPositiva] = 1
  return(out)
}

# Criterio 2 de proporcion de cobertura: Se utiliza un poligono convexo
.coverageIndex_2 = function(bitacoraAcustica, sp, grouped=TRUE) {
  # TO_DO: update isopData.old to be a parameter, inheriate it from estimarBiomasa
  sp = getSpeciesInfo(sp)$acustica
  data           = bitacoraAcustica[, c("lon", "lat", sp)]
  colnames(data) = tolower(colnames(data))
  
  data$area = .asignarIsoparalitorales(data)
  data      = data[complete.cases(data), ]
  
#   indexArea = data[, sp] > 0
  indexArea = tapply(data[, sp], data$area, sum)
  
  coverIndex = rep(0, len=length(unique(data$area)))
  names(coverIndex) = unique(data$area)

  for(i in names(indexArea)) {
    thisArea = as.numeric(i)
    tempData = as.matrix(data[data$area == thisArea, c("lon", "lat")])
    
    tempPolygon = chull(tempData)
    
    tempPolygon = tempData[c(tempPolygon, tempPolygon[1]), c("lon", "lat")]
    coords      = tempPolygon
    
    tempPolygon = SpatialPolygons(list(Polygons(list(Polygon(coords)), ID = 1)), 
                                  proj4string = CRS(proj4string(isopData.old)))
    
    tempArea = tempPolygon@polygons[[1]]@area
    
    isopIndex = match(thisArea, as.numeric(as.character(isopData.old$CODIGO)))
    isopArea  = isopData.old[isopIndex, 1]@polygons[[1]]@area
    
    coverIndex[i] = tempArea/isopArea
    
#     par(mar = c(2, 2, 1, 1))
#     thisCode <- match(thisArea, as.character(data3@data$CODIGO))
#     plot(data3@polygons[[thisCode]]@Polygons[[1]]@coords, type="l")
#     plot(tempPolygon, add = TRUE, border = "red")
#     points(tempData, pch = 19, cex = 0.5)
#     mtext(text = thisArea, side = 3, adj = 0.99, font = 2)
#     mtext(text = paste0("cov = ", round(tempArea/isopArea, 3)), side = 3, adj = 0.99, font = 2, line = -2)
  }
  
  out = coverIndex
  if(isTRUE(grouped)) {
    out = cut(out, breaks = c(-Inf, 0, 0.25, 0.5, 1.0), labels=c(0, 0.5,0.75,1))
    out = as.numeric(as.character(out))
  }

  names(out) = unique(data$area)
  
  return(out)
}

# Criterio 3 de proporcion de cobertura: Se utiliza un poligono convexo
.coverageIndex_3 = function(bitacoraAcustica, sp, grouped=FALSE){
  return(.coverageIndex_2(bitacoraAcustica, sp, grouped))
}

# Criterio 4 de proporcion de cobertura: Se utiliza un poligono concavo
.coverageIndex_4 = function(bitacoraAcustica, sp, grouped=TRUE) {
  # TO_DO: update isopData.old to be a parameter, inheriate it from estimarBiomasa
  sp = getSpeciesInfo(sp)$acustica
  data           = bitacoraAcustica[, c("lon", "lat", sp)]
  colnames(data) = tolower(colnames(data))
  
  data$area = .asignarIsoparalitorales(data)
  data      = data[complete.cases(data), ]
  
  #   indexArea = data[, sp] > 0
  indexArea = tapply(data[, sp], data$area, sum)
  
  coverIndex = rep(0, len=length(unique(data$area)))
  names(coverIndex) = unique(data$area)
  
  for(i in names(indexArea)[1:10]) {
    thisArea = as.numeric(i)
    tempData = as.matrix(data[data$area == thisArea, c("lon", "lat")])
    
    # Remove duplicated rows
    tempData = tempData[!duplicated(tempData),]
    if(is.vector(tempData))
      tempData = matrix(tempData, ncol = 2) else 
        tempData = as.matrix(tempData)
    
    if(!is.vector(tempData) && nrow(tempData) > 3){
      tempPolygon = .getArcs(ahull(x = tempData[,1] + runif(nrow(tempData), 0, 0.005), 
                                   y = tempData[,2] + runif(nrow(tempData), 0, 0.005), 
                                   alpha = 1))
    }else{
      tempPolygon = chull(tempData)
      
      tempPolygon = tempData[c(tempPolygon, tempPolygon[1]),]
      tempPolygon = list(x = tempPolygon[,1], y = tempPolygon[,2])
    }
    
    tempPolygon = data.frame(lon = tempPolygon$x,
                             lat = tempPolygon$y,
                             stringsAsFactors = FALSE)
    
    coords      = tempPolygon
    
    tempPolygon = SpatialPolygons(list(Polygons(list(Polygon(coords)), ID = 1)), 
                                  proj4string = CRS(proj4string(isopData.old)))
    
    tempArea = tempPolygon@polygons[[1]]@area

    isopIndex = match(thisArea, as.numeric(as.character(isopData.old$CODIGO)))
    isopArea  = isopData.old[isopIndex, 1]@polygons[[1]]@area
    
    coverIndex[i] = tempArea/isopArea
    
#     par(mar = c(2, 2, 1, 1))
#     thisCode <- match(thisArea, as.character(data3@data$CODIGO))
#     plot(data3@polygons[[thisCode]]@Polygons[[1]]@coords, type="l")
#     plot(tempPolygon, add = TRUE, border = "red")
#     points(tempData, pch = 19, cex = 0.5)
#     mtext(text = thisArea, side = 3, adj = 0.99, font = 2)
#     mtext(text = paste0("cov = ", round(tempArea/isopArea, 3)), side = 3, adj = 0.99, font = 2, line = -2)
  }
  
  out = coverIndex
  if(isTRUE(grouped)) {
    out = cut(out, breaks = c(-Inf, 0, 0.25, 0.5, 1.0), labels=c(0, 0.5,0.75,1))
    out = as.numeric(as.character(out))
  }
  
  names(out) = unique(data$area)
  
  return(out)
}

# Criterio 5 de proporcion de cobertura: Se utiliza un poligono convexo
.coverageIndex_5 = function(bitacoraAcustica, sp, grouped=FALSE){
  return(.coverageIndex_4(bitacoraAcustica, sp, grouped))
}

.getArcs <- function(x){
  onlyArc <- function (c, r, v, theta, ...) 
  {
    angles <- anglesArc(v, theta)
    seqang <- seq(angles[1], angles[2], length = 100)
    x <- c[1] + r * cos(seqang)
    y <- c[2] + r * sin(seqang)
    
    return(matrix(c(x, y), ncol = 2))
  }
  
  allArcs <- NULL
  arcs <- which(x$arcs[, 3] > 0)
  if (length(arcs) > 0) {
    for (i in arcs) {
      tempArcs <- onlyArc(x$arcs[i, 1:2], x$arcs[i, 3], x$arcs[i, 4:5], 
                          x$arcs[i, 6])
      
      allArcs <- rbind(allArcs, tempArcs)
    }
  }
  
  return(list(x = allArcs[,1], y = allArcs[,2]))
}

# calculate TS
.calculateTS = function(sp) {
  
  specie = getSpeciesInfo(sp) # get species information
  ts = getTS(sp)
  marcas = .createMarks(specie)
  thr = ts$thr
  b20 = ts$b20
  m   = ts$m
  
  if(length(thr)!=(length(b20)-1)) stop("Wrong number of b20 values.")
  thr = c(0, thr, Inf)
  xF = cut(marcas, breaks=thr, labels = FALSE)
  b20 = b20[xF]
  TS = m*log10(marcas) + b20
  return(TS)
}


.getDistLengthMap = function(data){
  
  datJuv = subset(data, subset = data$l <= 11.5)
  datAd1 = subset(data, subset = data$l > 11.5 & data$l <= 13.5)
  datAd2 = subset(data, subset = data$l > 13.5)
  
  labelsX=parse(text=paste(abs(c(84, 80, 76, 72)), "^o ", "*W", sep=""))
  labelsY=parse(text=paste(c(0, 5, 10, 15, 20), "^o ", "*S", sep=""))
  
  par(xaxs = "i", yaxs = "i")
  plot(datJuv$lon, datJuv$lat, pch = 19, type = "p", cex = 1, 
       col = "red", xlim = c(-84, -70), ylim = c(-20, 0), ylab = "",
       xlab = "", axes = FALSE)
  points(datAd1$lon, datAd1$lat, pch = 22, cex = 1, col = "blue", lwd = 2)
  points(datAd2$lon, datAd2$lat, pch = 3, cex = 0.8, col = "black", lwd = 2)
  map("worldHires", add=T, fill=T, col = "khaki1")
  axis(1, at = c(-84, -80, -76, -72), labels = 
         labelsX, cex.axis = 0.8, padj = -0.8)
  axis(2, at = c(0, -5, -10, -15, -20), labels = 
         labelsY, las = 2, cex.axis = 0.8, hadj = 0.7)
  mtext(text = "Latitud", side = 1, line = 2)
  mtext(text = "Latitud", side = 2, line = 2.5)
  legend(x = -75, y = -7, legend = c("Juvenil", "Adulta (12-13.5 cm)","Adulta (>13.5 cm)"), 
         pch = c(19, 22, 3), col = c("red", "blue", "black"), bty = "n")
  box()
  
  return(invisible())
  
}