R/CrunchIntersections.R

In mjkerrison/OzMap: Geographical Data for Australian Statistics

#' Calculate the Intersections
#'
#' An agglomeration of everything I did to find the intersections between the
#' LGAs and POAs.
#'
#' @param BigLGAobject The high-res LGA SPDF (see \code{\link{BigLGA}})
#' @param BigPOAobject The high-res POA SPDF (see \code{\link{BigPOA}})
#'
#' @export
CrunchIntersections <- function(BigLGAobject, BigPOAobject){


  # Utilities -----------------------------------------------------------------
  # For when we loop over the POAs/LGAs


  # To get some sp functions to work, we have to specify our coordinate system
  CoordRefSystem <- raster::crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  # I got this from somewhere on the internet, likely StackExchange, but can't recall where.
  # Note also that I don't fully understand the mathematical ins and outs of various specifications,
  # so you may take exception to this one.


  check.bounds <- function(range.A, range.B){
    # Is either of the (horizontal/vertical) extremities of A between the extremities of B?
    (range.B[1] <= range.A[1] & range.A[1] <= range.B[2]) | (range.B[1] <= range.A[2] & range.A[2] <= range.B[2])
    # If all four extremities of the bounding box of polygon A are outside that of polygon B,
    # then there's no chance they intersect.
  }


  crunch.intersection <- function(upper.poly, lower.poly){

    # Convert the polygons (LGA, POA) to SpatialPolygons objects so raster fns work
    upper.poly <- sp::SpatialPolygons(list(upper.poly), as.integer(1), CoordRefSystem)
    lower.poly <- sp::SpatialPolygons(list(lower.poly), as.integer(1), CoordRefSystem)

    # Note that the structue of the SpatialPolygonsDataFrame (in terms of the polygons) is:
    #   SPDF
    #       @polygons
    #                [[index]]*
    #                         @Polygons
    #                                  [[index]] - Formal class Polygon (various attr including coords)

    # SpatialPolygons wants the /collection/ of Polygons at an index of the SPDF (marked w/ *)

    # It's a laborious structure that took quite a while to figure out.

    # Quick check to see if we can skip or if we should call the slower(?) raster::intersect
    if(!any(

      # Check the y-bounds
      check.bounds(
        range(upper.poly@polygons[[1]]@Polygons[[1]]@coords[,1]),
        range(lower.poly@polygons[[1]]@Polygons[[1]]@coords[,1])
      ),

      # Check the x-bounds
      check.bounds(
        range(upper.poly@polygons[[1]]@Polygons[[1]]@coords[,2]),
        range(lower.poly@polygons[[1]]@Polygons[[1]]@coords[,2])
      )

    )){

      return(0)

    }

    # Otherwise, calculate the intersection between the two polygons as calculated by raster
    # (This returns a polygon itself)
    output <- raster::intersect(upper.poly, lower.poly)

    # Now beware: although the /bounds/ intersect, there's no guarantee that the actual /shapes/ do:
    if (is.null(output)) {

      return(0)

    } else {

      # Now calculate the area of this polygon: because we've got a SpatialPolygons object,
      # and our CRS specifies we're using long/lat, this is in square meters
      output <- raster::area(output)

      # Divide that by 1,000,000 to arrive at square kilometres, which is what comes with the shapefiles:

      return(output / 1000000)

      # And we return the ABSOLUTE VALUE so that, combined with those shapefiles, we can calculate
      # the percentage overlap either way.

    }

  }


  # Mapping over polygons -----------------------------------------------------

  # Use a progress bar from dplyr - this took ~40 mins on my computer
  progressbar.LGA <- dplyr::progress_estimated(nrow(BigLGAobject)*nrow(BigPOAobject))

  intersection.matrix <-

    # Map over the LGA polygons, returning rows,
    purrr::map_dfr(BigLGAobject@polygons, function(lga.poly){

      # which are constructed of columns corresponding to each POA polygon.
      purrr::map_dfc(BigPOAobject@polygons, function(poa.poly){

        # Tick the progress bar,
        progressbar.LGA$tick()$print()

        # Then return the km^2 intersection of this row (LGA) and column (POA)
        return(crunch.intersection(lga.poly, poa.poly))

      })

    })

  # Convert to matrix for size/efficiency (?)
  intersection.matrix <- as.matrix(intersection.matrix)

  # Name appropriately
  rownames(intersection.matrix) <- as.character(BigLGAobject$LGA_NAME16)
  colnames(intersection.matrix) <- as.character(BigPOAobject$POA_NAME16)

  # Voila
  return(intersection.matrix)

}