R/fun_vms_db_clean.R

In vmsbase: GUI Tools to Process, Analyze and Plot Fisheries Data

#' Neural Network Plotting - Internal Function
#'
#'
#' The \code{plotNet} is the internal function that implements the
#'  Neural Network Plotting routine.
#'
#' This function,  with a Neural Network, a Prediction Index and a Confusion Matrix, generated by \code{\link{gui_vms_met_pred}},
#'  plots the graphical overall result of the neural network prediction.
#'
#'
#' @return This function does not return a value.
#'
#' @param net The Neural Network
#' @param Dg The Prediction Index
#' @param ConfMat The Confusion Matrix
#'
#' @usage plotNet(net, Dg, ConfMat)
#'
#'
#' @references free text reference Pointers to the literature related to this object.
#' @seealso \code{\link{gui_vms_met_pred}}
#' @name plotNet

plotNet <- function(net, Dg, ConfMat) {
  TrErr <- net$Merror[, 1]
  VaErr <- net$Merror[, 2]

  min <- min(ConfMat)
  max <- max(ConfMat)
  yLabels <- rownames(ConfMat)
  xLabels <- colnames(ConfMat)
  title <- c("Confusion Matrix")
  # Red and green range from 0 to 1 while Blue ranges from 1 to 0
  ColorRamp <- rgb(
    seq(0, 1, length = 256), # Red
    seq(0, 1, length = 256), # Green
    seq(1, 0, length = 256)
  ) # Blue
  ColorLevels <- seq(min, max, length = length(ColorRamp))
  # Reverse Y axis
  reverse <- nrow(ConfMat):1
  yLabels <- yLabels[reverse]
  ConfMat <- ConfMat[reverse, ]


  # Define layout
  layout(matrix(data = c(1, 2, 1, 3), nrow = 2, ncol = 2), widths = c(4, 1, 4, 1), heights = c(4, 4))
  # Plot Error
  par(las = 2)
  plot(TrErr,
    type = "l",
    ylim = c(0, 1.1 * max(net$Merror)),
    lwd = 2, col = 2,
    xlab = "Iteration",
    ylab = "Mean Error"
  )
  lines(VaErr, lwd = 2, col = 4)
  legend(0, # 0.8*length(TrErr),
    y = 0.6 * max(net$Merror),
    legend = c("Training", "Validation"),
    fill = c(2, 4), col = c(2, 4), cex = 0.8
  )
  # Plot Conf Mat
  # Data Map
  par(mar = c(7, 11, 2, 4), las = 2)
  image(1:length(xLabels), 1:length(yLabels), t(ConfMat),
    col = ColorRamp, xlab = "",
    ylab = "", axes = FALSE, zlim = c(min, max)
  )
  if (!is.null(title)) {
    title(main = title)
  }
  axis(BELOW <- 1, at = 1:length(xLabels), labels = xLabels, cex.axis = 0.7)
  axis(LEFT <- 2,
    at = 1:length(yLabels), labels = yLabels, las = HORIZONTAL <- 1,
    cex.axis = 0.7
  )

  # Color Scale
  par(mar = c(3, 2.5, 2.5, 2))
  image(1, ColorLevels,
    matrix(data = ColorLevels, ncol = length(ColorLevels), nrow = 1),
    col = ColorRamp,
    xlab = "", ylab = "",
    xaxt = "n"
  )
}

StandardizeByCol <- function(xdata) {
  ncolx <- ncol(xdata)
  for (ij in 1:ncolx) {
    minx <- min(xdata[, ij], na.rm = TRUE)
    maxx <- max(xdata[, ij], na.rm = TRUE)
    xdata[, ij] <- (xdata[, ij] - minx) / (maxx - minx)
  }
  return(xdata)
}


#' Shape File Data Joining - Internal Function
#'
#'
#' The \code{Join2shp} is the internal function that implements the
#'  Shape File Data Joining routine.
#'
#' This function,  with a Grid Shape File, an Effort Count Vector and a File Destination Directory,
#'  (see \code{\link{gui_out_grid}}),
#'  creates a new Grid Shape File annotated with the Effort Count Vector.
#'
#'
#' @return This function does not return a value.
#'
#' @param shpfile The Original Grid Shape File
#' @param datavector The Effort Count Vector
#' @param dirdest The File Destination Directory
#'
#' @usage Join2shp(shpfile, datavector, dirdest)
#'
#'
#' @references free text reference Pointers to the literature related to this object.
#' @seealso \code{\link{gui_out_grid}}
#' @name Join2shp

Join2shp <- function(shpfile, datavector, dirdest) {
  sys_typ <- ifelse(.Platform$OS.type == "windows", "\\\\", "/")
  diror <- paste(unlist(strsplit(shpfile, sys_typ))[1:length(unlist(strsplit(shpfile, sys_typ))) - 1], collapse = sys_typ)
  if (diror == dirdest) {
    dirdest <- paste(dirdest, sys_typ, "grid", sep = "")
    dir.create(dirdest)
  }
  file_anm <- unlist(strsplit(unlist(strsplit(shpfile, sys_typ))[length(unlist(strsplit(shpfile, sys_typ)))], "[.]"))[1]

  otherfiles <- dir(diror)
  otherfiles <- otherfiles[grep(file_anm, otherfiles)]

  new_name <- ginput("Enter a name for\n the new shp file", title = "Shp File Name")

  for (iff in 1:length(otherfiles))
  {
    file.copy(
      from = paste(diror, sys_typ, otherfiles[iff], sep = "", collapse = ""),
      to = paste(dirdest, sys_typ, new_name, ".", unlist(strsplit(otherfiles[iff], "[.]"))[length(unlist(strsplit(otherfiles[iff], "[.]")))], sep = "", collapse = "")
    )
  }
  dbfdata <- read.dbf(paste(diror, sys_typ, file_anm, ".dbf", sep = ""), as.is = TRUE)

  dbfdata$Count <- datavector

  write.dbf(dbfdata, paste(dirdest, sys_typ, new_name, ".dbf", sep = "", collapse = ""))
}

# Original spline basis
H0 <- function(x) return(2 * x^3 - 3 * x^2 + 1)
H1 <- function(x) return(-2 * x^3 + 3 * x^2)
H2 <- function(x) return(x^3 - 2 * x^2 + x)
H3 <- function(x) return(x^3 - x^2)
# Derivatives
dH0 <- function(x) return(6 * x^2 - 6 * x)
dH1 <- function(x) return(-6 * x^2 + 6 * x)
dH2 <- function(x) return(3 * x^2 - 4 * x + 1)
dH3 <- function(x) return(3 * x^2 - 2 * x)

# by Josh O'Brien on http://stackoverflow.com/questions/9186496/determining-utm-zone-to-convert-from-longitude-latitude
long2UTM <- function(long) {
  (floor((long + 180) / 6) %% 60) + 1
}

# # Original spline basis
# H0 = function(x) return((2*x^3)-(3*x^2)+1)
# H1 = function(x) return(-(2*x^3)+(3*x^2)
# H2 = function(x) return((x^3)-(2*x^2)+x)
# H3 = function(x) return((x^3)-(x^2))
#
# # Derivatives
# dH0 = function(x) return((6*x^2) - (6*x))
# dH1 = function(x) return(-(6*x^2) + (6*x))
# dH2 = function(x) return((3*x^2) - (4*x) + 1)
# dH3 = function(x) return((3*x^2) - (2*x))


recu_dep <- function(xmin, xmax, ymin, ymax, resolut = 2, the_db = "") {
  if (the_db != "") {
    co_ce <- fn$sqldf("select count(*) from intrp where LON > `xmin` and LON < `xmax` and LAT > `ymin` and LAT < `ymax`", dbname = the_db)

    if (co_ce[1, 1] == 0) {
      cat("\n   -     Skipped block: x(", xmin, ",", xmax, ")*y(", ymin, ",", ymax, ")     -\n", sep = "")
    } else {
      xrange <- c(xmin, xmax)
      yrange <- c(ymin, ymax)
      # 	map("worldHires", xlim = xrange, ylim = yrange, bg = "darkorange2", col = "black", fill = T)
      # 	map.axes()
      l_x <- xmax - xmin
      l_y <- ymax - ymin
      if (l_x <= 0.25 | l_y <= 0.25) {
        cat("\n   -     New valid block: x(", xmin, ",", xmax, ")*y(", ymin, ",", ymax, ")     -\n", sep = "")
        deppoi <- fn$sqldf("select ROWID, * from intrp where LON > `xmin` and LON < `xmax` and LAT > `ymin` and LAT < `ymax`", dbname = the_db)

        cat("\n   -     Analyzing ", nrow(deppoi), " points     -\n\n", sep = "")

        bat_blo <- getNOAA.bathy(xmin - 0.1,
          xmax + 0.1,
          ymin - 0.1,
          ymax + 0.1,
          resolution = resolut
        )
        plot(bat_blo, image = T)
        points(deppoi[, "LON"], deppoi[, "LAT"], pch = 20, col = "firebrick")

        xlon <- rep(as.numeric(rownames(bat_blo)), length(as.numeric(colnames(bat_blo))))
        ylat <- rep(as.numeric(colnames(bat_blo)), each = length(as.numeric(rownames(bat_blo))))
        zdep <- as.numeric(bat_blo)
        cat("\n   -     Calculating Spline...     -", sep = "")
        SplineD <- o_Tps(cbind(xlon, ylat), zdep, lon.lat = TRUE)
        rm(bat_blo, zdep, xlon, ylat)
        cat("\n   -     Predicting depth", sep = "")
        if (nrow(deppoi) <= 10000) {
          cat(" ", sep = "")
          dept <- as.numeric(predict(SplineD, deppoi[, c("LON", "LAT")]))
          dep_v <- as.data.frame(cbind(deppoi[, "rowid"], deppoi[, "I_NCEE"], dept))
          sqldf("insert into p_depth select * from `dep_v`", dbname = the_db)
          rm(dept, dep_v)
          gc()
        } else {
          nPin <- ceiling(nrow(deppoi) / 10000)
          for (pi in 1:nPin)
          {
            cat(".", sep = "")
            r1 <- 10000 * (pi - 1) + 1
            r2 <- min(nrow(deppoi), r1 + 10000 - 1)
            dept <- as.numeric(predict(SplineD, deppoi[r1:r2, c("LON", "LAT")]))
            dep_v <- as.data.frame(cbind(deppoi[r1:r2, "rowid"], deppoi[r1:r2, "I_NCEE"], dept))
            sqldf("insert into p_depth select * from `dep_v`", dbname = the_db)
            rm(dept, dep_v)
            gc()
          }
        }
        cat(" - Completed!     -\n", sep = "")
        rm(SplineD)
        gc()
      } else {
        xmin_2 <- xmin + ((xmax - xmin) / 2)
        ymin_2 <- ymin + ((ymax - ymin) / 2)
        cat("\n   -     Splitting block...     -\n", sep = "")
        recu_dep(xmin, xmin_2, ymin, ymin_2, resolut = resolut, the_db)
        recu_dep(xmin_2, xmax, ymin, ymin_2, resolut = resolut, the_db)
        recu_dep(xmin, xmin_2, ymin_2, ymax, resolut = resolut, the_db)
        recu_dep(xmin_2, xmax, ymin_2, ymax, resolut = resolut, the_db)
      }
    }
  } else {
    cat("\n\n   -     Error! Bad Database File      -\n", sep = "")
  }
}

recu_dep_RDS <- function(bat_all, xmin, xmax, ymin, ymax, the_db = "", max_siz = 0.5, the_bbo) {
  if (the_db != "") {
    co_ce <- fn$sqldf("select count(*) from intrp where LON >= `xmin` and LON <= `xmax` and LAT >= `ymin` and LAT <= `ymax`", dbname = the_db)

    if (co_ce[1, 1] == 0) {
      cat("\n  -  Skipped block: x(", xmin, ",", xmax, ")*y(", ymin, ",", ymax, ")  -\n", sep = "")
    } else {
      xrange <- c(xmin, xmax)
      yrange <- c(ymin, ymax)
      l_x <- xmax - xmin
      l_y <- ymax - ymin
      if (l_x <= max_siz | l_y <= max_siz) {
        cat("\n  -  New valid block", sep = "")

        deppoi <- fn$sqldf("select ROWID, * from intrp where LON >= `xmin` and LON <= `xmax` and LAT >= `ymin` and LAT <= `ymax`", dbname = the_db)

        cat(" with ", nrow(deppoi), " points", sep = "")

        cur_rows <- which(as.numeric(rownames(bat_all)) >= xmin - 0.1 & as.numeric(rownames(bat_all)) <= xmax + 0.1)
        cur_cols <- which(as.numeric(colnames(bat_all)) >= ymin - 0.1 & as.numeric(colnames(bat_all)) <= ymax + 0.1)

        if (length(cur_rows) > 1 & length(cur_cols) > 1) {

          # cat("\n rows: ", length(cur_rows), " - cols:", length(cur_cols), "\n", sep = "")

          bat_blo <- bat_all[cur_rows, cur_cols]

          xlon <- rep(as.numeric(rownames(bat_blo)), length(as.numeric(colnames(bat_blo))))
          ylat <- rep(as.numeric(colnames(bat_blo)), each = length(as.numeric(rownames(bat_blo))))
          zdep <- as.numeric(bat_blo)
          #
          #           cat("\n xlon: ", xlon, "\n\n xlat: ", ylat, "\n\n zdep: ", zdep, "\n\n", sep = " ")
          #
          ##########
          blues <- c("lightsteelblue4", "lightsteelblue3", "lightsteelblue2", "lightsteelblue1")
          ba_bl <- as.bathy(cbind(xlon, ylat, zdep))
          par(mar = c(1, 1, 1, 1))
          plot(ba_bl, image = TRUE, land = TRUE, lwd = 0.1, bpal = list(c(0, max(zdep), "grey"), c(min(zdep), 0, blues)))
          plot(ba_bl, deep = 0, shallow = 0, step = 0, lwd = 0.4, add = TRUE)
          ##########

          #         plot(as.bathy(cbind(xlon, ylat, zdep)), image = T)

          points(deppoi[, "LON"], deppoi[, "LAT"], pch = 20, cex = 0.6, col = "firebrick")

          cat("  -  Calculating Spline...", sep = "")
          SplineD <- o_Tps(cbind(xlon, ylat), zdep, lon.lat = TRUE)
          #         cat(class(SplineD), "\n")
          rm(bat_blo, zdep, xlon, ylat)

          cat(" Predicting depth...", sep = "")
          if (nrow(deppoi) <= 10000) {
            dept <- as.numeric(predict(SplineD, deppoi[, c("LON", "LAT")]))
            dep_v <- as.data.frame(cbind(deppoi[, "rowid"], deppoi[, "I_NCEE"], dept))
            sqldf("insert into p_depth select * from `dep_v`", dbname = the_db)
            rm(dept, dep_v)
            gc()
          } else {
            nPin <- ceiling(nrow(deppoi) / 10000)
            for (pi in 1:nPin)
            {
              cat(".", sep = "")
              r1 <- 10000 * (pi - 1) + 1
              r2 <- min(nrow(deppoi), r1 + 10000 - 1)
              dept <- as.numeric(predict(SplineD, deppoi[r1:r2, c("LON", "LAT")]))
              dep_v <- as.data.frame(cbind(deppoi[r1:r2, "rowid"], deppoi[r1:r2, "I_NCEE"], dept))
              sqldf("insert into p_depth select * from `dep_v`", dbname = the_db)
              rm(dept, dep_v)
              gc()
            }
          }
          cat(" - Completed!  -\n", sep = "")
          map("world", xlim = the_bbo[2:1], ylim = the_bbo[4:3], col = "honeydew3", bg = "lightsteelblue1", fill = T)
          map.axes()

          xmin_2 <- xmin + ((xmax - xmin) / 2)
          abline(v = xmin_2, col = "firebrick")
          ymin_2 <- ymin + ((ymax - ymin) / 2)
          abline(h = ymin_2, col = "firebrick")

          rect(xmin, ymin, xmin_2, ymin_2, border = "firebrick")
          rect(xmin_2, ymin, xmax, ymin_2, border = "firebrick")
          rect(xmin, ymin_2, xmin_2, ymax, border = "firebrick")
          rect(xmin_2, ymin_2, xmax, ymax, border = "firebrick")

          rm(SplineD)
          gc()
        } else {
          cat(" outside downloaded bathymetry  -\n", sep = "")
        }
      } else {
        map("world", xlim = the_bbo[2:1], ylim = the_bbo[4:3], col = "honeydew3", bg = "lightsteelblue1", fill = T)
        map.axes()

        title("\n  -  Splitting Block...")

        xmin_2 <- xmin + ((xmax - xmin) / 2)
        abline(v = xmin_2, col = "firebrick")
        ymin_2 <- ymin + ((ymax - ymin) / 2)
        abline(h = ymin_2, col = "firebrick")

        rect(xmin, ymin, xmin_2, ymin_2, border = "firebrick")
        rect(xmin_2, ymin, xmax, ymin_2, border = "firebrick")
        rect(xmin, ymin_2, xmin_2, ymax, border = "firebrick")
        rect(xmin_2, ymin_2, xmax, ymax, border = "firebrick")

        recu_dep_RDS(bat_all, xmin, xmin_2, ymin, ymin_2, the_db, max_siz, the_bbo)
        recu_dep_RDS(bat_all, xmin_2, xmax, ymin, ymin_2, the_db, max_siz, the_bbo)
        recu_dep_RDS(bat_all, xmin, xmin_2, ymin_2, ymax, the_db, max_siz, the_bbo)
        recu_dep_RDS(bat_all, xmin_2, xmax, ymin_2, ymax, the_db, max_siz, the_bbo)
      }
    }
  } else {
    cat("\n\n   -     Error! Bad Database File      -\n", sep = "")
  }
}

#' Assign Area - Internal Function
#'
#'
#' The \code{Assign_Area} is the internal function that implements the
#'  Assign Area routine.
#'
#' This function,  with a VMS DB Track and a Sea Areas Shape File,
#'  (see \code{\link{gui_vms_db_are}}),
#'  finds the Sea Area that contain the VMS DB Track.
#'
#'
#' @return This function does not return a value.
#'
#' @param evnt The VMS Track Data
#' @param box The Sea Areas Shape File Box
#'
#' @usage Assign_Area(evnt, box)
#'
#'
#' @references free text reference Pointers to the literature related to this object.
#' @seealso \code{\link{gui_vms_db_are}}
#' @name Assign_Area

Assign_Area <- function(evnt, box) {
  Area <- max(findPolys(as.EventData(data.frame(EID = 1, X = mean(evnt[, 1]), Y = mean(evnt[, 2])), projection = "LL"), box)[, 2])
  if (length(Area) == 0) Area <- NA
  return(Area)
}



merge_source <- function(data_source1, data_source2, rnd_level, minover) {

  # Inspect sources
  uS1 <- unique(data_source1[, 1])
  nuS1 <- length(uS1)
  uS2 <- unique(data_source2[, 1])
  nuS2 <- length(uS2)

  #   #Set parameters
  #   npmin <- 30
  #   minover <- 30
  #   rnd_level <- 3

  # Select the sources with the highest number of vessels
  if (nuS1 < nuS2) {
    main_source <- data_source1
    secondary_source <- data_source2
  } else {
    main_source <- data_source2
    secondary_source <- data_source1
  }

  # Compute basic statistics
  uS1 <- unique(main_source[, 1])
  nuS1 <- length(uS1)
  uS2 <- unique(secondary_source[, 1])
  nuS2 <- length(uS2)

  matd <- matrix(data = Inf, nrow = nuS1, ncol = nuS2)
  rownames(matd) <- uS1
  colnames(matd) <- uS2
  iter <- 0

  for (i in 1:nuS1) {
    block_s1 <- main_source[which(main_source[, 1] == uS1[i]), ]
    # Spatial Overlap
    spat_data_source <- secondary_source[which((secondary_source$LON <= max(block_s1$LON)) &
      (secondary_source$LON >= min(block_s1$LON)) &
      (secondary_source$LAT <= max(block_s1$LAT)) &
      (secondary_source$LAT >= min(block_s1$LAT))), ]
    candidate_ID <- unique(spat_data_source$I_NCEE)
    nuC <- length(candidate_ID)

    for (j in 1:nuC) {
      block_s2 <- spat_data_source[which(spat_data_source[, 1] == candidate_ID[j]), ]

      # Semisincronicity
      overtimes_ss <- intersect(
        round(block_s1$DATE, rnd_level),
        round(block_s2$DATE, rnd_level)
      )
      if (length(overtimes_ss) > minover) {
        block_s1 <- block_s1[order(round(block_s1$DATE, rnd_level)), ]
        block_s2 <- block_s2[order(round(block_s2$DATE, rnd_level)), ]
        sub_block_s1 <- block_s1[which(round(block_s1$DATE, rnd_level) %in% overtimes_ss), ]
        sub_block_s2 <- block_s2[which(round(block_s2$DATE, rnd_level) %in% overtimes_ss), ]
        if (length(which(duplicated(round(sub_block_s1$DATE, rnd_level)) == TRUE)) > 0) {
          sub_block_s1 <- sub_block_s1[-which(duplicated(round(sub_block_s1$DATE, rnd_level)) == TRUE), ]
        }
        if (length(which(duplicated(round(sub_block_s2$DATE, rnd_level)) == TRUE)) > 0) {
          sub_block_s2 <- sub_block_s2[-which(duplicated(round(sub_block_s2$DATE, rnd_level)) == TRUE), ]
        }
        inter_blocks_dists <- gmt::geodist(sub_block_s1[, c("LAT")],
          sub_block_s1[, c("LON")],
          sub_block_s2[, c("LAT")],
          sub_block_s2[, c("LON")],
          units = "km"
        )
        matd[which(uS1 == uS1[i]), which(uS2 == candidate_ID[j])] <- median(inter_blocks_dists)
      }
    }
    iter <- iter + 1
    cat(iter, " of ", nuS1, " ", nuC, " candidates found", "\n")
  }
  return(matd)
}



o_Tps <- function(x, Y, m = NULL, p = NULL, scale.type = "range",
                  lon.lat = FALSE, miles = TRUE, ...) {
  x <- as.matrix(x)
  d <- ncol(x)
  if (is.null(p)) {
    if (is.null(m)) {
      m <- max(c(2, ceiling(d / 2 + 0.1)))
    }
    p <- (2 * m - d)
    if (p <= 0) {
      stop(" m is too small  you must have 2*m -d >0")
    }
  }
  Tpscall <- match.call()
  if (!lon.lat) {
    Tpscall$cov.function <- "Thin plate spline radial basis functions (Rad.cov) "
    Krig(x, Y,
      cov.function = Rad.cov, m = m, scale.type = scale.type,
      p = p, GCV = TRUE, ...
    )
  }
  else {
    # use different coding of the radial basis fucntions to use with great circle distance.
    Tpscall$cov.function <- "Thin plate spline radial basis functions (RadialBasis.cov) using great circle distance "
    Krig(x, Y,
      cov.function = stationary.cov, m = m, scale.type = scale.type,
      GCV = TRUE, cov.args = list(
        Covariance = "RadialBasis",
        M = m, dimension = 2, Distance = "rdist.earth",
        Dist.args = list(miles = miles)
      ), ...
    )
  }
}



#' Effort Count - Internal Function
#'
#'
#' The \code{CountMap} is the internal function that implements the
#'  Effort Count routine.
#'
#' This function,  with a VMS DB Track and a Sea Areas Shape File,
#'  (see \code{\link{gui_out_grid}}),
#'  computes the Effort Count Vector.
#'
#'
#' @return This function does not return a value.
#'
#' @param xy The VMS Fishing Point Data
#' @param GridPS The Sea Area Grid File
#'
#' @usage CountMap(xy, GridPS)
#'
#'
#' @references free text reference Pointers to the literature related to this object.
#' @seealso \code{\link{gui_out_grid}}
#' @name CountMap

CountMap <- function(xy, GridPS) {
  PP <- xy[, c(1:2)]
  GridCount <- matrix(0, length(unique(GridPS$PID)), 2)
  colnames(GridCount) <- c("cellID", "Count")
  if (nrow(PP) <= 100000) {
    Poi <- matrix(0, nrow(PP), 3)
    colnames(Poi) <- c("EID", "X", "Y")
    Poi[, 1] <- c(1:nrow(Poi))
    Poi[, 2] <- PP[, 1]
    Poi[, 3] <- PP[, 2]
    idPolys <- findPolys(Poi, GridPS)
    idTable <- table(idPolys[, 2])
    GridCount[, 1] <- 1:length(unique(GridPS$PID))
    GridCount[as.numeric(names(idTable)), 2] <- as.numeric(idTable)
  } else {
    nPP <- ceiling(nrow(PP) / 100000)
    #     GridCount = matrix(0,length(unique(GridPS$PID)),2)
    for (pi in 1:nPP) {
      cat(".", sep = "")
      r1 <- 100000 * (pi - 1) + 1
      r2 <- min(nrow(PP), r1 + 100000 - 1)
      Poi <- matrix(0, r2 - r1 + 1, 3)
      colnames(Poi) <- c("EID", "X", "Y")
      Poi[, 1] <- c(1:nrow(Poi))
      Poi[, 2] <- PP[r1:r2, 1]
      Poi[, 3] <- PP[r1:r2, 2]
      # GridPS = SpatialPolygons2PolySet(GSA)
      idPolys <- findPolys(Poi, GridPS) # Supporta fino a 100.000 punti
      idTable <- table(idPolys[, 2])
      GridCount_part <- matrix(0, length(unique(GridPS$PID)), 2)
      colnames(GridCount_part) <- c("cellID", "Count")
      GridCount_part[, 1] <- 1:length(unique(GridPS$PID))
      GridCount_part[as.numeric(names(idTable)), 2] <- as.numeric(idTable)
      GridCount[, 2] <- GridCount[, 2] + GridCount_part[, 2]
    }
  }
  return(GridCount[, 2])
}