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