R/pop_dasymap.R

In ecological-cities/home2park: Spatial Provision of Urban Parks

#'Perform dasymetric mapping of (rasterised) population data
#'
#'Distribute population data per census block across inhabitable land as defined by `land_relative_density`.
#'Raster data provided should be a snapshot of the population for a single time period (e.g. a year),
#'created with other package functions (e.g. `rasterise_pop()`, `rasterise_buildings()`).
#'If there are multiple years present in the intermediate datasets, the user will have to extract the relevant list element for
#'the single year of interest (see examples). All input data should have a similar projected coordinate reference system
#'specific to the target area.
#'
#'@param pop_polygons `sf` polygons of the population data for a single time period.
#'May be extracted from the `pop_polygons` list element generated by `rasterise_pop()`.
#'@param pop_perblock_count Raster of population count per census block for a single time period.
#'May be extracted from the `pop_count_rasters` list element generated by `rasterise_pop()`.
#'@param pop_perblock_density Raster of population density per census block for a single time period.
#'May be extracted from the `pop_density_rasters` list element generated by `rasterise_pop()`.
#'@param land_relative_density Raster of the relative density (e.g. suitability, habitability) of land for a single time period.
#'May be extracted from results generated by `rasterise_buildings()` or `rasterise_landuse()`.
#'@param filename character (optional). Export output raster to disk.
#'@param overwrite logical. Argument passed to `terra::writeRaster()`. If `TRUE`, `filename` is overwritten.
#'@param ... Other arguments passed to `terra::writeRaster()`.
#'
#'@return Raster of population density, based on relative density values of the land defined in `land_relative_density`.
#'
#'@import sf
#'@import checkmate
#'@importFrom dplyr bind_cols select filter
#'@importFrom glue glue
#'@importFrom terra vect rast rasterize
#'@importFrom rlang .data
#'
#'@examples
#' \dontrun{
#' data(pop_sgp) # population census block polygons
#' data(landuse_sgp) # land use polygons
#'
#'
#' # transform to projected crs
#' pop_sgp <- sf::st_transform(pop_sgp, sf::st_crs(32648))
#' landuse_sgp <- sf::st_transform(landuse_sgp, sf::st_crs(32648))
#'
#'
#' # merge all census blocks for chosen year (2020) into single multi-polygon
#' # function requires that polygons are merged
#' city_boundaries <- pop_sgp %>%
#'    dplyr::filter(year == 2020) %>%
#'    sf::st_union() %>%
#'    sf::st_as_sf() %>%
#'    smoothr::fill_holes(threshold = units::set_units(1, 'km^2'))  %>% # clean up
#'    smoothr::drop_crumbs(threshold = units::set_units(1, 'km^2'))  %>%
#'    sf::st_make_valid()
#'
#' buildings <- get_buildings_osm(place = city_boundaries,
#'                                date = as.Date('2021-01-01')) %>%
#'    mutate(year = 2020)
#'
#'
#' # rasterise population, landuse & buildings
#' pop_rasters <- rasterise_pop(pop_sgp,
#'                              census_block = "subzone_n",
#'                              pop_count = "pop_count")
#'
#' landuse_rasters <- rasterise_landuse(landuse_sgp,
#'                                      land_use = 'lu_desc',
#'                                      subset = c('1' = 'RESIDENTIAL',
#'                                                 '2' = 'COMMERCIAL & RESIDENTIAL',
#'                                                 '3' = 'RESIDENTIAL WITH COMMERCIAL AT 1ST STOREY',
#'                                                 '4' = 'RESIDENTIAL / INSTITUTION'),
#'                                      sf_pop = pop_sgp,
#'                                      match_landuse_pop = 'recent')
#'
#' buildings_rasters <- rasterise_buildings(buildings,
#'                                          proxy_pop_density = 'levels',
#'                                          year = 'year',
#'                                          sf_pop = pop_sgp,
#'                                          sf_landuse = landuse_sgp,
#'                                          match_buildings_pop = 'closest')
#'
#'
#' # finally, perform dasymetric mapping on selected year (2020)
#' popdens_raster <- pop_dasymap(pop_polygons = pop_rasters$pop_polygons[[2]],
#'                               pop_perblock_count = pop_rasters$pop_count[[2]],
#'                               pop_perblock_density = pop_rasters$pop_density[[2]],
#'                               land_relative_density = buildings_rasters[[2]],
#'                               filename = 'buildings_popdensity.tif',
#'                               wopt = list(gdal=c('COMPRESS=LZW')))
#' }
#'
#'@export
pop_dasymap <- function(pop_polygons, pop_perblock_count, pop_perblock_density, land_relative_density,
    filename = NULL, overwrite = TRUE, ...) {


    # Error checking ------------------

    coll <- checkmate::makeAssertCollection()

    # file formats for polygons
    checkmate::assertTRUE(!st_is_longlat(pop_polygons) & !is.null(st_crs(pop_polygons)), add = coll)  # must be projected crs
    checkmate::assertTRUE(all(st_is_valid(pop_polygons)), add = coll)  # all features must be valid

    # all crs similar
    checkmate::assertTRUE(all(st_crs(pop_polygons) == st_crs(pop_perblock_count) & st_crs(pop_polygons) ==
        st_crs(pop_perblock_density) & st_crs(pop_polygons) == st_crs(land_relative_density)))

    # file paths
    checkmate::assert_character(filename, min.len = 1, any.missing = FALSE, all.missing = FALSE,
        null.ok = TRUE, add = coll)

    checkmate::reportAssertions(coll)


    # Calculations ------------------

    # 1) Get relative densities per census block polygon multiply uniform pop density per
    # census block with relative density grid (suitable pixels & no. of levels for pop data to
    # spread across)
    pop_perblock_rel_density <- pop_perblock_density * land_relative_density


    # 2) Zonal sum of pop_perblock_rel_density (per census block) - later serves as
    # non-negative weights tt sum to 1 within each
    pop_perblock_rel_density_sum <- terra::extract(pop_perblock_rel_density, vect(pop_polygons),
        fun = "sum", na.rm = TRUE)

    colnames(pop_perblock_rel_density_sum) <- c("ID", "sum")

    # overwrite, assume same order
    pop_perblock_rel_density_sum <- pop_polygons %>%
        dplyr::bind_cols(pop_perblock_rel_density_sum) %>%
        dplyr::select(-.data$ID)


    # rasterize
    pop_perblock_rel_density_sum <- terra::rasterize(vect(pop_perblock_rel_density_sum), land_relative_density,
        field = "sum")


    # multiply pop_perblock_rel_density & pop_perblock_count: to adjust weights used to make
    # the sum of values equal to the total count in each polygon divide by
    # pop_perblock_rel_density_sum: set of non-negative weights tt sum to unity within each
    # polygon prod(res(pop_perblock_rel_density_sum)) # no need to divide by area per pixel
    results <- pop_perblock_rel_density * pop_perblock_count/pop_perblock_rel_density_sum


    # export
    if (!is.null(filename)) {

        terra::writeRaster(results, filename = filename, overwrite = overwrite, ...)
    }

    terra::tmpFiles(remove = TRUE)
    gc()

    return(results)
}