R/rasterise_landuse.R

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

#'Rasterise residential land use polygons
#'
#'Rasterise land use zones (`sf_landuse`) with reference to population data (`sf_pop`), if supplied.
#'Multiple (a list of) output rasters will be generated if land use data for multiple years are present.
#'If the land use data includes unnecessary land use classes, the argument `subset` and `landuse`
#'allows the user to first subset the data to polygons defined as 'residential' land use.
#'
#'If population data (`sf_pop`) is supplied, its raster(s) must have been previously generated in the `tempdir()` using
#'the function `rasterise_pop()`. Each year in the population data is matched to a specific year in the land use data
#'(using the helper function `matchyear()`). Land use raster(s) will then
#'be associated with a specific population census year, and removed if there are no matching population census years.
#'Rasterised land use zones will be masked away (i.e., convert pixels to `NA`)
#'at areas with no (zero) population data (for the respective matching year).
#'
#'@param sf_landuse `sf` polygons of the land use zones.
#'Data should be in a projected coordinate reference system.
#'@param subset Named vector containing the residential land use classes of interest to retain in column `land_use`.
#'Each element is named according to the integer that will be used
#'to represent this class in the output raster. Defaults to `NULL` (integers will be alphabetically
#'assigned to land use classes).
#'@param land_use character. Specify column name of the land use classes.
#'@param year character. Specify column name for the year within `sf_landuse`
#'and `sf_pop` (if provided). Defaults to `'year'`. Column data should be numeric.
#'@param sf_pop (optional) `sf` polygons containing the population census data with column containing the census year.
#'If absent, the output (land use rasters) will not be associated with specific population census year(s).
#'@param match_landuse_pop character. Type of matching between `sf_landuse` and `sf_pop`, passed on to `match` argument in function `matchyear()`.
#'Either `'exact'`, `'closest'`, `'recent'` or `'soonest'`. Defaults to `'recent'`, i.e.
#'assumes that land use is pre-planned (and followed by population changes), rather than a result of post-monitoring (e.g. remotely-sensed).
#'@param dir_rastertemplate character. Filepath to the raster used to define the pixel resolution, extent, nrow, ncol of
#'the output  raster; object is passed to the '`y`' argument in `terra::rasterize()`.
#'Defaults to the template raster generated by the function `rasterise_pop()` within `dir_processing` (see next argument).
#'@param dir_processing character. Directory to get intermediate files generated in previous steps (e.g. rasterised population data).
#'Defaults to `tempdir()`.
#'@param dir_export character. File path to directory to export output raster(s) to. Defaults to `tempdir()`.
#'Set to `NULL` if you do not wish to export output for subsequent processing.
#'@param overwrite logical. Argument passed to `terra::writeRaster()`. Defaults to `TRUE`.
#'@param wopt list. Argument passed to `terra::writeRaster()`.
#'@param ... Other arguments passed to `terra::writeRaster()`.
#'
#'@return List of raster files. Zero values are converted to `NA`.
#'Each raster is also an intermediate file exported to the directory as defined by `dir_export`
#'(defaults to `tempdir()`) for further processing.
#'
#'@import sf
#'@import checkmate
#'@importFrom glue glue
#'@importFrom dplyr filter mutate
#'@importFrom terra vect rast rasterize
#'@importFrom rlang .data
#'
#'@examples
#' \dontrun{
#' data(pop_sgp)
#' data(landuse_sgp)
#'
#'
#' # 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))
#'
#'
#' # run function
#' 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'),
#'                                      year = 'year',
#'                                      sf_pop = pop_sgp,
#'                                      match_landuse_pop = 'recent')
#' }
#'
#'@export
rasterise_landuse <- function(sf_landuse, subset = NULL, land_use = NULL, year = "year", sf_pop = NULL,
    match_landuse_pop = "recent", dir_rastertemplate = NULL, dir_processing = tempdir(), dir_export = tempdir(),
    overwrite = TRUE, wopt = list(gdal = c("COMPRESS=LZW")), ...) {

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

    coll <- checkmate::makeAssertCollection()

    # file formats
    checkmate::assertTRUE(!st_is_longlat(sf_landuse) & !is.null(st_crs(sf_landuse)), add = coll)  # must be projected crs
    # checkmate::assertTRUE(all(st_is_valid(sf_landuse)), add = coll) # all features must be
    # valid (too time consuming)

    # colnames
    checkmate::assert_subset(land_use, choices = colnames(sf_landuse), empty.ok = FALSE, add = coll)
    checkmate::assert_subset(year, choices = colnames(sf_landuse), empty.ok = TRUE, add = coll)

    # land use classes
    checkmate::assert_character(subset, min.len = 1, any.missing = FALSE, all.missing = FALSE,
        null.ok = TRUE, add = coll)
    checkmate::assert_subset(subset, choices = unique(sf_landuse[[land_use]]), empty.ok = TRUE,
        add = coll)  # as vector

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

    checkmate::reportAssertions(coll)


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

    # IF POP DATA IS PROVIDED
    if (!is.null(sf_pop)) {

        # match pop data to most recent landuse data
        sf_pop <- matchyear(data = sf_pop, data_tomatch = sf_landuse, year = year, match = match_landuse_pop)

        # remove landuse years not required (no matching pop years)
        sf_landuse <- sf_landuse %>%
            dplyr::filter(.data[[year]] %in% unique(sf_pop$year_match))

        if (nrow(sf_landuse) == 0) {
            # none fulfill criteria!
            stop("Error: No land use data fulfills the criteria specified in 'match_landuse_pop'!")
        }

    }


    # SUBSET LAND USE
    if (!is.null(subset)) {

        sf_landuse <- sf_landuse %>%
            dplyr::filter(.data[[land_use]] %in% subset)
    }

    # create new col with integers representing land use classes
    sf_landuse <- sf_landuse %>%
        dplyr::mutate(residential = match(.data[[land_use]], subset))  # if subset = NULL, just col w NAs


    # GET RASTER TEMPLATE IF USER SUPPLIED
    if (!is.null(dir_rastertemplate)) {
        raster_template <- rast(dir_rastertemplate)
    } else {
        raster_template <- rast(glue::glue("{dir_processing}/popdensity_raster-template.tif"))
    }


    # RASTERIZE
    results <- list()  # create empty list for output

    # If there is a year column in sf_landuse
    if (!is.null(unique(sf_landuse[[year]])))
        {

            years_landuse <- unique(sf_landuse[[year]])

            # each landuse year
            for (i in seq_along(years_landuse)) {

                sf_landuse_subset <- sf_landuse %>%
                  dplyr::filter(.data[[year]] == years_landuse[i])

                # rasterize & export to temp directory
                landuse_raster <- terra::rasterize(terra::vect(sf_landuse_subset), raster_template,
                  field = "residential")

                landuse_raster[landuse_raster == 0] <- NA  # convert 0s to NAs


                # IF POP DATA IS PROVIDED, remove areas not covered by pop census blocks (for matching
                # census yr)
                if (!is.null(sf_pop)) {

                  years_popmatch <- unique(sf_pop$year[sf_pop$year_match == years_landuse[i]])
                  years_popmatch <- na.omit(years_popmatch)

                  results[[i]] <- list()

                  for (j in seq_along(years_popmatch)) {

                    pop_raster_matching <- rast(glue::glue("{dir_processing}/popdensity-by-census-block_{years_popmatch[j]}.tif"))

                    # mask away landuse zones not within relevant pop subzones
                    results[[i]][[j]] <- terra::mask(landuse_raster, pop_raster_matching)
                    names(results[[i]][[j]]) <- glue::glue("landuse-residential-{years_landuse[i]}_for-pop-{years_popmatch[j]}")

                    # export
                    if (!is.null(dir_export)) {
                      terra::writeRaster(results[[i]][[j]], filename = file.path(glue::glue("{dir_export}/{names(results[[i]][[j]])}.tif")),
                        overwrite = overwrite, wopt = wopt, ...)
                    }
                    terra::tmpFiles(remove = TRUE)  # remove temp files to clear memory
                    rm(j)
                  }


                  # POP DATA NOT PROVIDED
                } else {

                  results[[i]] <- landuse_raster
                  names(results[[i]]) <- glue::glue("landuse-residential-{years_landuse[i]}")

                  # export
                  if (!is.null(dir_export)) {
                    terra::writeRaster(results[[i]], filename = file.path(glue::glue("{dir_export}/{names(results[[i]])}.tif")),
                      overwrite = overwrite, wopt = wopt, ...)
                  }
                  terra::tmpFiles(remove = TRUE)  # remove temp files to clear memory
                  rm(i)
                }

                rm(landuse_raster)
            }


            # If no year column in sf_landuse
        }  # else{
    # results[[1]] <- terra::rasterize(terra::vect(sf_landuse), raster_template, field =
    # 'residential') results[[1]][results[[1]] == 0] <- NA # convert 0s to NAs # export
    # names(results[[1]]) <- 'landuse-residential' terra::writeRaster(results[[1]], filename =
    # file.path(glue::glue('{dir_processing}/{names(results[[1]])}')), overwrite = TRUE, wopt
    # = list(gdal=c('COMPRESS=LZW'))) rm(i) }

    terra::tmpFiles(remove = TRUE)  # remove temp files to clear memory
    return(unlist(results))
}