# 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.