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