R/quadtree.R

In tomauer/sptrees:

#' Quadtree
#'
#' @export
quadtree <- function(spts,
                     buffer,
                     max_points_per_bin,
                     max_bin_area,
                     min_bin_area,
                     rotation = NA,
                     n_cores = 1) {

  # rotate the points before binning
  # but set 0 to NA
  if(!is.na(rotation)) {
    if(rotation == 0) {
      rotation <- NA
    }
  }

  if(!is.na(rotation)) {
    spts <- maptools::elide(spts,
                          rotate = rotation,
                          center = c(rgeos::gCentroid(spts)@coords[1],
                                     rgeos::gCentroid(spts)@coords[2]))
    raster::crs(spts) <- sp::CRS(sp::proj4string(spts))
  }

  # create stripped down object
  xy <- data.frame(x = spts@coords[, 1],
                   y = spts@coords[, 2])

  binDepths <- c(0)
  binParents <- c(0)
  binCorners <- list(c(min(xy$x) - buffer,
                       min(xy$y) - buffer,
                       max(xy$x) + buffer,
                       max(xy$y) + buffer))

  pointBins <- rep(1, nrow(xy))

  checkPoints <- function(bcs, xy) {
    cpFun <- function(i) {
      if(xy[i,]$x > bcs[1] &
         xy[i,]$y > bcs[2] &
         xy[i,]$x < bcs[3] &
         xy[i,]$y < bcs[4]) {
        return(TRUE)
      } else {
        return(FALSE)
      }
    }

    containsPoints <- parallel::mclapply(1:nrow(xy),
                                         cpFun,
                                         mc.preschedule = TRUE,
                                         mc.set.seed = TRUE,
                                         mc.cores = n_cores)

    return(unlist(containsPoints))
  }

  divide <- function(binNo) {

    # calculate new bin dimensions and store old bin dimensions
    newDiv <- (binCorners[[binNo]][1:2] + binCorners[[binNo]][3:4]) / 2
    oldDiv <- binCorners[[binNo]]

    # remove the old bin, since we're storing
    binCorners[[binNo]] <<- NULL

    lapply(1:4, function(i) {

      newBinNo <- length(binCorners) + 1

      if(i == 1) {
        binCorners[[newBinNo]] <<- c(oldDiv[1:2], newDiv)
      } else if(i == 2) {
        binCorners[[newBinNo]] <<- c(oldDiv[1], newDiv[2], newDiv[1], oldDiv[4])
      } else if(i == 3) {
        binCorners[[newBinNo]] <<- c(newDiv, oldDiv[3:4])
      } else if(i == 4) {
        binCorners[[newBinNo]] <<- c(newDiv[1], oldDiv[2:3], newDiv[2])
      }



      if(!(sum(pointBins == binNo) == 0)) {
        pts <- checkPoints(binCorners[[newBinNo]], xy[pointBins == binNo, ])
        pointBins[pointBins == binNo] <<- ifelse(pts, newBinNo, binNo)

        # square kilometers
        # is this method for calculating area safe?
        # calculated in square meters
        binArea <- ((binCorners[[newBinNo]][3] -
                       binCorners[[newBinNo]][1]) *
                      (binCorners[[newBinNo]][4] -
                         binCorners[[newBinNo]][2]))

        # points per sq km
        #binDensity <- sum(pts)/(binArea/1000000)

        # sum(pts) > 100
        # binDensity > 50
        # 1200 sq km would be a good lower

        #print(sum(pts))
        #print(binArea)

        # need to add a minimum bin area


        if( (sum(pts) > max_points_per_bin | binArea > max_bin_area) &
            !(binArea < (min_bin_area * 4))) {
          divide(newBinNo)
        }
      }
    })

    return(binCorners)
  }

  qt <- divide(1)

  qtdf <- do.call(rbind.data.frame, qt)
  names(qtdf) <- c("xmin", "ymin", "xmax", "ymax")
  qtdf$ids <- 1:nrow(qtdf)

  ID <- row.names(qtdf)

  square <- cbind(qtdf$xmin, qtdf$ymax, qtdf$xmax, qtdf$ymax, qtdf$xmax,
                  qtdf$ymin, qtdf$xmin, qtdf$ymin, qtdf$xmin, qtdf$ymax)

  polys <- SpatialPolygons(mapply(function(poly, id) {
    xy <- matrix(poly, ncol=2, byrow=TRUE)
    Polygons(list(Polygon(xy)), ID=id)
  }, split(square, row(square)), ID),
  proj4string = sp::CRS(sp::proj4string(spts)))

  tdspolydf <- SpatialPolygonsDataFrame(polys, qtdf)
  tdspolydf$area <- rgeos::gArea(tdspolydf, byid = T)

  return(tdspolydf)
}