R/path-auxiliar.R

In imarpe/path: Pelagic Assessment Tools for Hydroacoustics

# get TS from data
getTS = function(sp) { 
  
  if(!is.character(sp)) sp = deparse(substitute(sp))
  if(!(sp %in% names(.TSdata))) stop(sprintf("unknown species %s.", sp))
  
  output = .TSdata[[sp]]
  
  class(output) = c("path.ts", class(output))
  
  return(output)
}

# filter a sub-range of Latitudes
.filterLatitude = function(base, lat) {
  if(!is.null(lat)) {
    ind = !is.na(cut(base$lat, breaks=lat))
    base = base[ind,  ]
    return(base)
  } else return(base) 
}


# Genera un vector con valores relativos a la suma de los valores del vector
.relative = function(x)
  return(x/sum(x, na.rm = TRUE))

# Extrae valores de latitud a partir de uno de area isoparalitoral
.extraeLatitude = function(x) {
  latitude = as.numeric(substring(text = x, first = nchar(x) - 2, last = nchar(x) - 1))
  return(latitude)
}

# Calcula distancias entre dos puntos
.distancia = function(x1, x2, y1, y2, unidad = "mn"){
  
  lon1 = -x1*60*cos(-y1*pi/180)
  lon2 = -x2*60*cos(-y2*pi/180)
  lat1 = -y1*60
  lat2 = -y2*60
  
  out = sqrt((lon1 - lon2)^2 + (lat1 - lat2)^2)
  
  if(unidad == "mn")
    out = out*1852
  
  return(out)
}


asignarRegion = function(data) {
  region = character(nrow(data))
  region[data$lat>-16] = "nc"
  region[data$lat<=-16] = "sur"
  return(region)
}



#Funcion para asignar areas isoparalitorales
.asignarIsoparalitorales = function(data, colLon = NA, colLat = NA, old = TRUE) {
  
  if(is.na(colLon))
    colLon = match("lon", tolower(colnames(data)))
  
  if(is.na(colLat))
    colLat = match("lat", tolower(colnames(data)))
  
  data = data.frame(lon = data[,colLon], lat = data[,colLat], stringsAsFactors = FALSE)
  isoCodes = numeric(nrow(data))
  
  # Read points and convert to SpatialPointsDataFrame object
  points = data
  
  isoCodes[!complete.cases(points)] = NA
  coords = points[!is.na(isoCodes),]
  coordinates(coords) = ~ lon + lat
  
  # Read polygon of 'Areas isoparalitorales' and convert to SpatialPolygonsDataFrame object 
  peru.isoparalit = if(old) isopData.old else isopData.new
  
  proj4string(coords) = sp::CRS(proj4string(peru.isoparalit))
  
  isoCodes[!is.na(isoCodes)] = as.numeric(as.matrix(sp::over(coords, peru.isoparalit)))
  
  return(isoCodes)
}


#Asigna los valores de a y b si es que no hay muestreo biológico o biométrico
.assignWeightParameters = function(baseCalas, sp) {
  # limits from north to south
  # to be deleted
  if(is.null(baseCalas[["a", exact=TRUE]])) baseCalas$a = getSpeciesInfo(sp)$a  
  if(is.null(baseCalas[["b", exact=TRUE]])) baseCalas$b = getSpeciesInfo(sp)$b  
  return(baseCalas)
}

# Función para convertir la base de datos cruda de crucero al formato clásico
.convertirCruceroMF = function (base, sp,  ...) {
  
  
  #Transformando Coordenadas
  base$lon_gi = abs(trunc(base$lon))
  base$lon_mi = abs(round((base$lon + base$lon_gi)*60, 2))
  base$lat_gi = abs(trunc(base$lat))
  base$lat_mi = abs(round((base$lat + base$lat_gi)*60, 2))
  
  #Datos de cada lance
  header1 = aggregate(captura ~ lance + buque + sp, FUN = unique, data=base, na.action=na.pass)
  header1$captura = as.numeric(sapply(header1$captura, mean, na.rm = TRUE))
  header1 = aggregate(captura ~ lance + buque, FUN=sum, data=header1, na.action=na.pass)
  header1$w_cap_lance_kg = header1$captura
  header1$captura = NULL
  
  #Filtra por especie
  base = base[base$sp == sp, ]
  
  header = aggregate(cbind(fecha, lat, lon, lon_gi, lon_mi, lat_gi, lat_mi, dc,
                           w_muestra_g, w_cap_anch_kg) ~ lance + buque, 
                     data = base, FUN = head, na.action=na.pass, n = 1)
  header = merge(header, header1, sort=FALSE)
  
  # Convertir header2 a modo numerico, a excepcion de la columna "fecha" y "buque"
  for(i in which(!(colnames(header) %in% c("fecha", "buque", "lance"))))
    header[, i] = suppressWarnings(as.numeric(as.character(header[, i])))
  
  #Frecuencia de tallas muestreadas

  calaBuque = paste(base$buque, base$lance, sep="-")
  
  freq = tapply(base$frec, INDEX = list(calaBuque, base$l), sum)
  freq[is.na(freq)] = 0
  
  #Crea marcas totales segun especie
  specie = getSpeciesInfo(sp)
  marks = .createMarks(specie)
  
  #Completa frecuencia de tallas no muestreadas
  freqmarks = matrix(0, nrow = nrow(freq), ncol = length(marks))
  colnames(freqmarks) = marks
  freqmarks[, colnames(freq)] = freq

  n = rowSums(freqmarks)
  
  xbaseFreqs = do.call(rbind, strsplit(rownames(freq), "-"))
  
  baseFreqs = data.frame(buque = xbaseFreqs[,1], lance = xbaseFreqs[,2], 
                         n = n, freqmarks, check.names = FALSE)

  #Base Transformada y filtrada
  base = merge(header, baseFreqs, sort=FALSE)
  rownames(base) = seq_len(nrow(base))
  
  # Calcular proporción de cobertura por area isoparalitoral y asignar a una columna 'Prop'
  base$prop = base$w_cap_anch_kg/base$w_cap_lance_kg
  base$area = .asignarIsoparalitorales(base[, c("lon", "lat")])
  
  return(base = base)
  
}

.convertirCruceroED = function(base, sp, ...) {
    
  # Convertir nombres de filas de base a minusculas
  rownames(base) = tolower(rownames(base))
  # Convertir base a modo numerico, a excepcion de la columna "fecha"
  base = as.data.frame(t(base), stringsAsFactors = FALSE)
  for(i in which(!(colnames(base) %in% c("fecha", "buque"))))
    base[,i] = suppressWarnings(as.numeric(base[,i]))

  marcasCols = grep(x = colnames(base), pat = "^[0-9]")
  marks = as.numeric(colnames(base)[marcasCols])
  
  base2 = base[, marcasCols]
  base2[is.na(base2)] = 0
  frec = as.numeric(as.matrix(base2))
  l = rep(marks, each=nrow(base2))
  
  base2 = data.frame(lance=base$lance, buque=base$buque, l=l, frec=frec)
  
  base1 = base[, -marcasCols]
#   
#   #Crear marcas según especie
#   specie = getSpeciesInfo(sp)
#   marks = .createMarks(specie)
#   
  # Calcular las coordenadas de cada cala y adjuntarlas en columnas respectivas
  base1$lon = -(base1$lon_g + base1$lon_m/60)
  base1$lat = -(base1$lat_g + base1$lat_m/60)
  
  # Calcular proporción de cobertura??? por area isoparalitoral y asignar a una columna 'Prop'
  base1$prop = base1$w_cap_anch_kg/base1$w_cap_lance_kg
  base1$area = .asignarIsoparalitorales(base1[, c("lon", "lat")])
  
  base = merge(x = base1, y = base2, sort = FALSE)
  ind = base$frec == 0
  base = base[!ind, ]

  base$sp = "anchoveta"
  base$peso = NA # TO_DO: parse 'peso.medio' for ED

  base$peso.muestra = base$w_muestra_g 
  base$captura = base$w_cap_anch_kg
  base$dc = kali::getDistance(data = base, ref = coast_noIsland)/1.852
#   base$w_muestra_g[is.na(base$w_muestra_g)] = 0 # TO_DO: check data values (error)

  return(base)
}

.checkFile = function(directorio, file) {
  if(is.null(file)) return(NULL)
  filedoc   = file.path(directorio, file)
  if(!file.exists(filedoc))
    stop(sprintf("El archivo %s no existe.", filedoc))
  return(filedoc)
}


.applySampleSizeThreshold = function(baseCalas, umbralIndividuos = 150){
    
  filasValidas  = which(baseCalas$n >= umbralIndividuos)
  
  # Volver a unir en baseCalas
  baseCalas = baseCalas[filasValidas, ]
  
  # TO_DO: print calas and buques
  if(length(filasValidas)>0)
      cat("Descartando calas:", baseCalas$lance[which(!filasValidas)], "\n")  
  
  return(baseCalas)
    
}

getDistArea <- function(data, coordNames = c("lon_m", "lat_m"),
                        lonGrid_params = c(-85, -70, 0.25), latGrid_params = c(-20, 0, 0.25), 
						convFactor = 60^2, 
                        units = "mn2", projData = NULL, xlim = NULL, ylim = NULL,
                        plot = FALSE, ...){
  
  # Set griding function
  gridPoints <- function(data, lonGrid_params, latGrid_params, remove.duplicated = TRUE){
    # 'data': A data.frame with just lon-lat columns
    
    as.nc <- function(x, ...){
      return(as.numeric(as.character(x, ...)))
    }
    
    lonGrid <- seq(lonGrid_params[1], lonGrid_params[2], lonGrid_params[3])
    latGrid <- seq(latGrid_params[1], latGrid_params[2], latGrid_params[3])
    
    gridData <- expand.grid(lon = lonGrid, lat = latGrid)
    gridData <- data.frame(ID = 1:nrow(gridData), gridData)
    
    griddedLon <- cut(data[,1], 
                      breaks = c(lonGrid - lonGrid_params[3]/2, max(lonGrid) + lonGrid_params[3]/2),
                      labels = lonGrid)
    griddedLat <- cut(data[,2], 
                      breaks = c(latGrid - latGrid_params[3]/2, max(latGrid) + latGrid_params[3]/2),
                      labels = latGrid)
    
    griddedPoints <- data.frame(lon = as.nc(griddedLon), lat = as.nc(griddedLat), 
                                stringsAsFactors = FALSE)
    griddedPoints$ID <- gridData$ID[match(with(griddedPoints, paste0(lon, lat)),
                                          with(gridData, paste0(lon, lat)))]
    
    if(isTRUE(remove.duplicated)){
      griddedPoints <- griddedPoints[!duplicated(with(griddedPoints, paste0(lon, lat))),]
    }
    
    return(list(griddedPoints = griddedPoints, gridData = gridData, 
                ngrids = sum(!duplicated(with(griddedPoints, paste0(lon, lat))))))
  }
  
  # Griding points
  data_gridded <- gridPoints(data = data[,coordNames],
                                lonGrid_params = lonGrid_params, 
                                latGrid_params = latGrid_params, 
                                remove.duplicated = TRUE)
  
  # Get only the data.frame of points
  data <- data_gridded$griddedPoints
  
  # Skip if some of data have minus than 3 rows
  # if(nrow(data) < 3)
  #   next
  
  # Get distribution area
  distArea <- data_gridded$ngrids*lonGrid_params[3]*latGrid_params[3]*convFactor
  distArea <- round(distArea, 1)
  
  
  # Plot distribution area 
  data$z <- rep(1, nrow(data))
  
  # Convert data.frame to raster object
  data_raster <- data[,c("lon", "lat", "z")]
  
  coordinates(data_raster) <- ~ lon + lat
  
  projection(data_raster) <- projData
  
  gridded(data_raster) <- TRUE  
  
  data_raster <- raster(data_raster)
  
  # Plot raster
  if(isTRUE(plot)){
    if(is.null(xlim))
      xlim <- lonGrid_params[1:2]
    
    if(is.null(ylim))
      ylim <- latGrid_params[1:2]
    
    image(data_raster, xlim = xlim, ylim = ylim, ...)
  }
    
  
  return(list(distribution_area = distArea,
              units = units,
              data = data,
              data_raster = data_raster))
}



getDistArea2 <- function(data, grilla, ...){

latitudes  = seq(-19, -3-grilla, by=grilla)
longitudes = seq(-84, -69-grilla, by=grilla)

grupos.lat = list()
valores    = matrix(NA, ncol = length(longitudes), nrow = length(latitudes))

for(i in seq_along(latitudes)){
  grupos.lat[[i]]=subset(data, data$lat_m > latitudes[i] &
                           data$lat_m < latitudes[i] + grilla)
  for(j in seq_along(longitudes)){
    tempData = subset(grupos.lat[[i]], grupos.lat[[i]]$lon_m > longitudes[j] &
               grupos.lat[[i]]$lon_m < longitudes[j] + grilla)
    
    if(sum(tempData$anc, na.rm = TRUE) == 0) { valores[i,j] = NA }
      else { valores[i,j] = 1 }
  }
}

x1     = grilla * 60
areamn = x1 * x1

axisZ   = as.vector(valores)
outtemp = sum(axisZ, na.rm = TRUE)
out = outtemp * areamn
  
  return(out)

  }


.getIndices = function(data, stock, dist.costa, lat.max, lat.min, gridsize, ...){
  
  names(data) = tolower(names(data))
  
  if(!(c("anc") %in% names(data))){
  
	x1 = which(names(data) == "ancp")
	x2 = which(names(data) == "ancg")

	data$anc = rowSums(as.data.frame(data[ , c(x1, x2)]))

  }
  
  dataList = list()
  
  if(stock == "nc") dataList[[1]] = subset(data, subset = data$lat_m >= -16)
  if(stock == "sur") dataList[[1]] = subset(data, subset = data$lat_m < -16)
  if(stock == "total") {
	  dataList[[1]] = subset(data, subset = data$lat_m >= -16)
	  dataList[[2]] = subset(data, subset = data$lat_m < -16)
	  dataList[[3]] = data
  }

  indMatrix = matrix(NA, ncol = 9, nrow = length(dataList))
  
  for(i in seq_along(dataList)){
  
	  data1 = dataList[[i]]
	  
	  if(nrow(data1) == 0) {  indicadores = matrix(rep(NA, 9), nrow=1)
	  
	  } else {
	  
	  names(data1) = tolower(names(data1))
	  data1 = data1[!is.na(data1$lat_m), ]
	  data1 = data1[!is.na(data1$lon_m), ]
	  data1 = data1[!is.na(data1$anc), ]
	  posiciones = data1[,c("lon_m","lat_m")]
	  #posiciones = na.omit(posiciones)

	  #- Convert VMS data to SpatialPolygons
	  spTa              = SpatialPoints(data.frame(posiciones))
	  proj4string(spTa) = CRS("+proj=longlat")
	  spTa.proj         = spTransform(spTa, CRS("+proj=utm +zone=18 ellips=WGS84"))
	  
	  #- Read shapefile of Peru
	  Peru              = as(PER_ADM0, "SpatialPolygons")
	  proj4string(Peru) = CRS("+proj=longlat")
	  Peru.proj         = spTransform(Peru, CRS("+proj=utm +zone=18 ellips=WGS84"))
	  dists = gDistance(spgeom1 = spTa.proj, spgeom2=Peru.proj,byid=T) #
	  data1$DC 		  = t(dists*0.00053996) # convirtiendo de metros a millas nauticas
	  
	  data1 = subset(data1, subset = data1$lat_m < lat.max & data1$lat_m > lat.min)
	  data1 = subset(data1, subset = data1$DC <= dist.costa)
	  
	  Plus   = nrow(data1[data1$anc!=0,])
	  Null   = nrow(data1[data1$anc==0,])
	  ISO    = (Plus/(Plus+Null))*100
	  sAPlus = mean(data1$anc[data1$anc!=0], na.rm=T)
	  sA     = mean(data1$anc, na.rm=T)
	  	  
	  CG_dc      = sum(data1$DC*data1$anc, na.rm = TRUE)/sum(data1$anc, na.rm = TRUE)
	  Inertia_dc = (sum(((data1$DC - CG_dc)^2)*data1$anc, na.rm = TRUE))/sum(data1$anc, na.rm = TRUE)
	  
	  CG_lat      = sum(data1$lat_m*data1$anc, na.rm = TRUE)/sum(data1$anc, na.rm = TRUE)
	  Inertia_lat = (sum(((data1$lat_m - CG_lat)^2)*data1$anc, na.rm = TRUE))/sum(data1$anc, na.rm = TRUE)
	  
	  projData <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
	  
	  #distArea <- getDistArea(data = data1, coordNames = c("lon_m", "lat_m"), projData = projData)
	  distArea <- getDistArea2(data = data1, grilla = gridsize)
	  Gini = gini(data1$anc, rep(1,length(data1$anc)))
	  #Discon = DCdensity(data1$anc[data1$anc!=0], sAPlus)
	  
	  #indicadores = matrix(c(distArea, ISO, sA, sAPlus, CG_dc, Inertia_dc, CG_lat, Inertia_lat, Gini, Discon), nrow=1)
	  indicadores = matrix(c(distArea, ISO, sA, sAPlus, CG_dc, Inertia_dc, CG_lat, Inertia_lat, Gini), nrow=1)

	}
	
  indMatrix[i, ] = indicadores
  
  }
	
  colnames(indMatrix) <- c("Area_dist", "ISO", "SA","SAplus","CG_DC", "Inertia_DC", "CG_Lat", "Inertia_Lat", "Gini")
  #colnames(indMatrix) <- c("Area_dist", "ISO", "SA","SAplus","CG_DC", "Inertia_DC", "CG_Lat", "Inertia_Lat", "Gini", "Discontinuity")
  
  if(stock == "nc") rownames(indMatrix) <- c("NC")
  if(stock == "sur") rownames(indMatrix) <- c("Sur")
  if(stock == "total") rownames(indMatrix) <- c("NC", "Sur", "Total")

  return(indMatrix)
  
}