R/demand_noise_reg.R
In ecoservR/ecoserv_tool: Natural Capital Mapping Toolkit (UK)
Defines functions
demand_noise_reg
Documented in
demand_noise_reg
#########################################################################
### Noise regulation demand model ###
### Adapted from EcoServ ###
### 10-11-2021 ###
### Sandra Angers-Blondin and Karl Busdieker ###
#########################################################################
# NOTES
# This script uses the classified Base Map to produce a raster of the demand for noise regulation within a study area.
# See Ecoserv technical report for descriptions, assumptions and caveats.
# 2021 update: roads, airports and motorways are now built into the package
# https://www.nature.scot/snh-research-report-954-ecoserv-gis-v33-toolkit-mapping-ecosystem-services-gb-scale
##########################################################################################################################
#' Noise Regulation Demand Model
#'
#' Runs the noise regulation ecosystem service model, generating demand scores based on three indicators: distance to noise sources (airports, railways and major roads), population, and health. (Specific indicators can be omitted by setting the appropriate weight to 0.)
#' @param x A basemap, in a list of sf tiles or as one sf object. Must have attributes "HabCode_B", "housePop", "health".
#' @param studyArea The boundaries of the site, as one sf object. The final raster will be masked to this shape. For best results this shape should be smaller than the basemap (which should be buffered by typically 300 m - 1km to avoid edge effects).
#' @param res Desired resolution of the raster. Default is 5 m. Range recommended is 5-10m.
#' @param distances Distance threshold (in meters) for noise attenuation for the noise sources: motorways, dual carriageways, primary roads, railways, and airports. Must be a named vector of numeric values in meters; defaults to c("motorway" = 800,"dual" = 600,"primary" = 550,"railway" = 650,"airport" = 1500). Defra's Strategic Noise Mapping (https://www.gov.uk/government/publications/strategic-noise-mapping-2019) may help set sensible defaults.
#' @param local Radius (m) for focal statistics at local range. Default is 300 m.
#' @param indicators Logical; should raw indicators (before transformation into z-scores and rescaling) be saved to the project folder? Default TRUE.
#' @param indic_weights A numeric vector of length 3, with weights for distance to roads, population and health, respectively. Default to c(1, 0.5, 0.5) so that population and health combined has equal weight to noise sources.
#' @param projectLog The RDS project log file generated by the wizard app and containing all file paths to data inputs and model parameters
#' @param runtitle A customised title you can give a specific model run, which will be appended to your project title in the outputs. If comparing a basemap to an intervention map, we recommend using "pre" and "post", or a short description of the interventions, e.g. "baseline" vs "tree planting".
#' @param save Path to folder where outputs will be saved. By default a folder will be created using your chosen run title, prefixed by "services_". Do not use this argument unless you need to save the outputs somewhere else.
#' @return Two rasters with demand scores: one with raw scores (arbitrary units), and one rescaled 0-100 (where 100 is maximum demand for the area).
#' @export
#'
demand_noise_reg <- function(x = parent.frame()$mm,
studyArea = parent.frame()$studyArea,
res = 5,
distances = c(
"motorway" = 800,
"dual" = 600,
"primary" = 550,
"railway" = 650,
"airport" = 1500
),
local = 300,
indicators = TRUE,
indic_weights = c(1, 0.5, 0.5),
projectLog = parent.frame()$projectLog,
runtitle = parent.frame()$runtitle,
save = NULL){
## Setup ----
timeA <- Sys.time() # start time
# Create output directory automatically if doesn't already exist
if (is.null(save)){
save <- file.path(projectLog$projpath,
paste0("services_", runtitle))
if (!dir.exists(save)){
dir.create(save)
}
} else {
# if user specified their own save directory we check that it's ok
if(!dir.exists(save) | file.access(save, 2) != 0){
stop("Save directory doesn't exist, or you don't have permission to write to it.")}
}
# Create a temp directory for scratch files
scratch <- file.path(projectLog$projpath,
"ecoservR_scratch")
if(!dir.exists(scratch)){
dir.create(scratch)
}
## Initial checks ----
attributes <- c("HabCode_B", "housePop", "health")
# Is the indicator argument an expected value (TRUE or FALSE)?
if(!is.logical(indicators)) {
warning("indicators argument must be TRUE or FALSE; using TRUE by default")
indicators <- TRUE
}
# Set (and create if needed) the indicator folder
if (indicators) {
indicator_path <- file.path(save, "indicators") # set the indicator path
if (!dir.exists(indicator_path)){dir.create(indicator_path)}
}
# Warning and resulting to default if weights entered incorrectly
try(if(length(indic_weights) != 3 | !is.numeric(indic_weights)){
warning("Weights must be a numeric vector of length 3 for this model. Using the default (noise:1, population: 0.5, health: 0.5). Specify a weight of 0 if you want to ignore a dataset.", call. = FALSE)
indic_weights <- c(1, 0.5, 0.5)
}
)
# Are the weights specified correctly?
if (is.numeric(indic_weights)){
if (length(indic_weights) != 3){ # if only one weight specified we assign it to each layer
message("Weights should have 3 values for this model. Defaulting to weights of c(1, 0.5, 0.5).")
indic_weights <- c(1,1,1,1)
}
} else stop("Weights must be a numeric vector of length 3 for this model. Please refer to documentation. Specify a weight of 0 if you want to ignore a dataset.")
# Check the distance thresholds
if (!all(names(distances) == c("motorway", "dual", "primary", "railway", "airport")) | !is.numeric(distances)){
stop("Distance thresholds set incorrectly. Must be a named vector: make sure that names match documentation. Alternatively remove the argument to use default settings.")
}
# if mm is stored in list, combine all before proceeding
if (isTRUE(class(x) == "list")){
x <- do.call(rbind, x) %>% sf::st_as_sf()
x <- checkgeometry(x, "POLYGON")
# NOT using rbindlist here because only keeps the extent of the first tile
}
# Are any attributes missing?
if (any(!attributes %in% names(x))){
stop("Could not find attribute(s) in your basemap: ",
paste(attributes[which(!attributes %in% names(x))], collapse = ", "),
". Please add or rename the attributes.")
}
### Keep only fields we need
message("Preparing basemap and study area.")
x <- x[attributes]
# Create a buffer around the study area (to clip roads with some margin)
SAbuffer <- sf::st_buffer(sf::st_geometry(studyArea), 500) %>% sf::st_as_sf() %>% sf::st_transform(27700)
# Create raster with same properties as mastermap to use as a template
r <- raster::raster() # create empty raster
raster::crs(r) <- sp::CRS(SRS_string = "EPSG:27700") # hard-coding datum to preserve CRS
raster::extent(r) <- raster::extent(x) # set same extent as the basemap
raster::res(r) <- res # set resolution
################################# DISTANCE RASTER ##############################
if (indic_weights[1] > 0){ # only if user wants this indicator
# NOTE - Higher scores here mean an area is further from a noise source, so scores are inverted to mean that higher
# scores represent greater need for noise regulation
message("Importing noise sources (OS VectorMap) from EcoservR")
##### Importing OS VectorMap data #####
suppressWarnings({
airports <- readRDS(system.file("extdata/vectormap/airports.RDS", package = "ecoservR")) %>% sf::st_set_crs(27700)
railways <- readRDS(system.file("extdata/vectormap/railways.RDS", package = "ecoservR")) %>% sf::st_set_crs(27700)
roads <- readRDS(system.file("extdata/vectormap/roads.RDS", package = "ecoservR")) %>% sf::st_set_crs(27700)
})
#### Clip them to the buffered study area
airports <- sf::st_intersection(airports, SAbuffer)
railways <- sf::st_intersection(railways, SAbuffer)
roads <- sf::st_intersection(roads, SAbuffer)
##### Subsetting roads and buffering sources #####
# Create layers for noise sources
# (ONLY if there are occurrences in the study area, otherwise create empty raster)
message("Preparing VectorMap data for noise distance analysis")
motorways <- dplyr::filter(roads, classification == "Motorway") %>% # subset of roads that are motorways
sf::st_buffer(15) %>%
checkgeometry() # buffer and convert to polygon
dualcarriageways <- dplyr::filter(roads, classification == "A Road") %>% # subset of roads that are dual carriageways
sf::st_buffer(5) %>%
checkgeometry() # buffer and convert to polygon
roadsA <- dplyr::filter(roads, classification == "Minor Road") %>% # subset of A roads
sf::st_buffer(5) %>%
checkgeometry() # buffer and convert to polygon
airports <- airports %>%
sf::st_buffer(500) %>%
checkgeometry()
railways <- railways %>%
sf::st_buffer(5) %>%
checkgeometry()
##### 4 - Calculating noise source distances #####
message("Creating noise distance raster (indicator 1 of 3)")
# calculate distance to noise sources, writing resulting rasters to scratch folder:
# here the maxdist argument could become a reactive argument defined by thresholds from the defra noise dataset
dist_motor <- distance_from_source(motorways, r, "motorways", distances[["motorway"]])
dist_roadsA <- distance_from_source(roadsA, r, "roadsA", distances[["primary"]])
dist_dualcarriageways <- distance_from_source(dualcarriageways, r, "dualcarriageways", distances[["dual"]])
dist_rail <- distance_from_source(railways, r, "railways", distances[["railway"]])
dist_air <- distance_from_source(airports, r, "airports", distances[["airport"]])
### TO DO: REVISE DISTANCES SET ABOVE AND BELOW (THRESHOLDS FOR EACH SOURCE)
## This was actual distances but the score must be inversely proportional.
## From ecoservGIS:converted to log10 score and then to inverse score
# also added check because there might not be any features of a kind, and empty raster causes error in calculation
dist_motor <- if (raster::hasValues(dist_motor)){1 - (log10(dist_motor + 1) / log10(distances[["motorway"]]))} else {r}
dist_roadsA <- if (raster::hasValues(dist_roadsA)){1 - (log10(dist_roadsA + 1) / log10(distances[["primary"]]))} else {r}
dist_dualcarriageways <- if (raster::hasValues(dist_dualcarriageways)){1 - (log10(dist_dualcarriageways + 1) / log10(distances[["dual"]]))} else {r}
dist_rail <- if (raster::hasValues(dist_rail)){1 - (log10(dist_rail + 1) / log10(distances[["railway"]]))} else {r}
dist_air <- if (raster::hasValues(dist_air)){1 - (log10(dist_air + 1) / log10(distances[["airport"]]))} else {r}
##### 5 - Summing distance rasters #####
# Sum rasters:
# NB: have to use stack function (originally brick, caused a bug)
distances <- raster::stack(dist_motor, dist_roadsA, dist_dualcarriageways, dist_rail, dist_air) # stack rasters
rm(dist_motor, dist_roadsA, dist_dualcarriageways, dist_rail, dist_air)
if (terra::hasValues(distances)){
distancescore <- raster::calc(distances, fun = sum, na.rm = TRUE) #sum scores
# score could be slightly negative (precision), so forcing to zero
distancescore[distancescore < 0] <- 0
} else {
# if no noise sources at all in study area, set raster to 0
distancescore <- r
terra::values(distancescore) <- 0
}
# save raw indicator
if (indicators){
distancescore <- raster::writeRaster(distancescore,
filename = file.path(indicator_path,
paste(projectLog$title,runtitle, "noise_reg_dist_to_noise_indic.tif", sep="_")),
overwrite = TRUE)
message("Distance to noise sources indicator saved.")
}
rm(airports, motorways, railways, roads, roadsA, dualcarriageways)
} else {
distancescore <- r # if user wishes to omit this indicator we create empty raster
}
### POP AND/OR HEALTH
## If at least the pop and/or health indicators are selected (non-zero weights), create the houses
if (indic_weights[2] > 0 | indic_weights[3] > 0){
# Create the spatial points with socioeconomic data for each household
# We can't go straight from polygon to raster, because then a big house that spans two raster cells would
# have its population counted twice. By converting to points first, we make sure that each house is only represented once. If more than house house falls within a raster cell, their population gets summed.
houses <- dplyr::filter(x, HabCode_B == "J360") %>% # filter houses from basemap
checkgeometry("POLYGON") %>% # need to be single-part before converting to centro
sf::st_centroid() # get the centroid of each house
}
# Population raster -------------------------------------------------------
if (indic_weights[2] > 0){ # only if user wants this indicator
# Population size (sum) within local search distance
# NOTE - Higher population scores mean there is more societal need
message("INDICATOR 2 of 3: population density")
## Create the population raster by summing individual house populations in each cell
popscore <- raster::rasterize(houses, r, field = "housePop", fun = "sum", na.rm = TRUE,
filename = file.path(scratch, "airpurif_pop_temp"), overwrite = TRUE)
message("...calculating population in ", local, " m radius")
# Calculate the summarised local scores
popscore <- focalScore(popscore, radius = local, type = "sum")
# ## OLD FOCAL STATS
# # Create weight matrices based on the search radius for focal stats
# # (automatically considers the res of the raster to calculate distance)
# window_local <- raster::focalWeight(r, local, "circle") # search window for the longer dist
# window_local[window_local > 0] <- 1 # replacing weights by 1 (we want full sum, not mean)
#
#
# message("...calculating population in ", local, " m radius")
#
# # Calculate the summarised local scores
# popscore <- raster::focal(popscore, w = window_local, na.rm = TRUE, pad = TRUE)
# Save raw indicator if need be
if (indicators){
raster::writeRaster(popscore, file.path(indicator_path, # we define the path
paste(projectLog$title,runtitle, "noise_reg_pop_indic.tif", sep = "_")),
overwrite = TRUE)
}
} else {
popscore <- r
}
# Health indicator --------------------------------------------------------
if (indic_weights[3] > 0){ # only if user wants this indicator
# Health deprivation score (mean) within local search distance
# NOTE - Higher health scores mean an area is more deprived, i.e. there is more societal need for noise regulation
message("INDICATOR 3 of 3: health score")
# Create the health raster from the houses data we subsetted in previous indicator
healthscore <- raster::rasterize(houses, r, field = "health", fun = mean, na.rm = TRUE,
filename = file.path(scratch, "airpurif_health_temp"),
overwrite = TRUE)
rm(houses)
# Calculate the summarised local scores
healthscore <- focalScore(healthscore, radius = local, type = "mean")
### OLD FOCAL STATS
# healthscore <- raster::focal(healthscore, w = window_local,
# fun = function(x){mean(x[window_local != 0], na.rm=TRUE)
# # we select only the values in the circular window (discarding values where window is 0), and perform a mean calculation on this subset rather than the full rectangular extent of the matrix. This ensure the correct denominator while discarding NAs.
# },
# pad = TRUE) # pad edges with NAs so raster has full extent)
# Save raw indicator if need be
if (indicators){
raster::writeRaster(healthscore, file.path(indicator_path, # we define the path
paste(projectLog$title, runtitle, "noise_reg_health_indic.tif", sep = "_")),
overwrite = TRUE)
}
} else {
healthscore <- r
}
# RESCALE INDICATORS ------------------------------------------------------
# Each indicator gets transformed into a z-score (mean-centered and deviance-scaled), and then reprojected on a 1-100 scale
distancescore <- rescale_indicator(distancescore)
popscore <- rescale_indicator(popscore)
healthscore <- rescale_indicator(healthscore)
# COMBINE INDICATORS ------------------------------------------------------
# All indicators should now be between 0-1 and we can sum them.
message("Combining all indicators...")
final <- sum(
indic_weights[1]*distancescore,
indic_weights[2]*popscore,
indic_weights[3]*healthscore,
na.rm = TRUE # important otherwise only gives raster for where all 4 indicators have values
)
rm(distancescore, popscore, healthscore)
# Replace NA with 0
final[is.na(final)] <- 0
# Mask to study area
final <- raster::mask(final, studyArea)
message("Saving final scores.")
raster::writeRaster(final,
filename = file.path(save, paste(projectLog$title, runtitle,
"noise_reg_demand.tif",
sep="_")),
overwrite = TRUE)
## and a rescaled version
maxvalue <- max(raster::values(final), na.rm = TRUE)
raster::writeRaster(final/maxvalue *100,
filename = file.path(save, paste(projectLog$title, runtitle,
"noise_reg_demand_rescaled.tif", sep="_")),
overwrite = TRUE)
timeB <- Sys.time() # start time
# write performance to log
projectLog$performance[["dem_noise"]] <- as.numeric(difftime(
timeB, timeA, units="mins"
))
updateProjectLog(projectLog) # save revised log
# Delete all the stuff we don't need anymore
on.exit({
rm(r, final, maxvalue)
cleanUp(scratch)
message("Noise regulation demand model finished. Process took ", round(difftime(timeB, timeA, units = "mins"), digits = 1), " minutes. Please check output folder for your maps.")
})
return({
## returns the objects in the global environment
invisible({
projectLog <<- projectLog
})
})
}
ecoservR/ecoserv_tool documentation built on April 5, 2025, 1:49 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions demand_noise_reg
Documented in demand_noise_reg
#########################################################################
### Noise regulation demand model ###
### Adapted from EcoServ ###
### 10-11-2021 ###
### Sandra Angers-Blondin and Karl Busdieker ###
#########################################################################
# NOTES
# This script uses the classified Base Map to produce a raster of the demand for noise regulation within a study area.
# See Ecoserv technical report for descriptions, assumptions and caveats.
# 2021 update: roads, airports and motorways are now built into the package
# https://www.nature.scot/snh-research-report-954-ecoserv-gis-v33-toolkit-mapping-ecosystem-services-gb-scale
##########################################################################################################################
#' Noise Regulation Demand Model
#'
#' Runs the noise regulation ecosystem service model, generating demand scores based on three indicators: distance to noise sources (airports, railways and major roads), population, and health. (Specific indicators can be omitted by setting the appropriate weight to 0.)
#' @param x A basemap, in a list of sf tiles or as one sf object. Must have attributes "HabCode_B", "housePop", "health".
#' @param studyArea The boundaries of the site, as one sf object. The final raster will be masked to this shape. For best results this shape should be smaller than the basemap (which should be buffered by typically 300 m - 1km to avoid edge effects).
#' @param res Desired resolution of the raster. Default is 5 m. Range recommended is 5-10m.
#' @param distances Distance threshold (in meters) for noise attenuation for the noise sources: motorways, dual carriageways, primary roads, railways, and airports. Must be a named vector of numeric values in meters; defaults to c("motorway" = 800,"dual" = 600,"primary" = 550,"railway" = 650,"airport" = 1500). Defra's Strategic Noise Mapping (https://www.gov.uk/government/publications/strategic-noise-mapping-2019) may help set sensible defaults.
#' @param local Radius (m) for focal statistics at local range. Default is 300 m.
#' @param indicators Logical; should raw indicators (before transformation into z-scores and rescaling) be saved to the project folder? Default TRUE.
#' @param indic_weights A numeric vector of length 3, with weights for distance to roads, population and health, respectively. Default to c(1, 0.5, 0.5) so that population and health combined has equal weight to noise sources.
#' @param projectLog The RDS project log file generated by the wizard app and containing all file paths to data inputs and model parameters
#' @param runtitle A customised title you can give a specific model run, which will be appended to your project title in the outputs. If comparing a basemap to an intervention map, we recommend using "pre" and "post", or a short description of the interventions, e.g. "baseline" vs "tree planting".
#' @param save Path to folder where outputs will be saved. By default a folder will be created using your chosen run title, prefixed by "services_". Do not use this argument unless you need to save the outputs somewhere else.
#' @return Two rasters with demand scores: one with raw scores (arbitrary units), and one rescaled 0-100 (where 100 is maximum demand for the area).
#' @export
#'
demand_noise_reg <- function(x = parent.frame()$mm,
studyArea = parent.frame()$studyArea,
res = 5,
distances = c(
"motorway" = 800,
"dual" = 600,
"primary" = 550,
"railway" = 650,
"airport" = 1500
),
local = 300,
indicators = TRUE,
indic_weights = c(1, 0.5, 0.5),
projectLog = parent.frame()$projectLog,
runtitle = parent.frame()$runtitle,
save = NULL){
## Setup ----
timeA <- Sys.time() # start time
# Create output directory automatically if doesn't already exist
if (is.null(save)){
save <- file.path(projectLog$projpath,
paste0("services_", runtitle))
if (!dir.exists(save)){
dir.create(save)
}
} else {
# if user specified their own save directory we check that it's ok
if(!dir.exists(save) | file.access(save, 2) != 0){
stop("Save directory doesn't exist, or you don't have permission to write to it.")}
}
# Create a temp directory for scratch files
scratch <- file.path(projectLog$projpath,
"ecoservR_scratch")
if(!dir.exists(scratch)){
dir.create(scratch)
}
## Initial checks ----
attributes <- c("HabCode_B", "housePop", "health")
# Is the indicator argument an expected value (TRUE or FALSE)?
if(!is.logical(indicators)) {
warning("indicators argument must be TRUE or FALSE; using TRUE by default")
indicators <- TRUE
}
# Set (and create if needed) the indicator folder
if (indicators) {
indicator_path <- file.path(save, "indicators") # set the indicator path
if (!dir.exists(indicator_path)){dir.create(indicator_path)}
}
# Warning and resulting to default if weights entered incorrectly
try(if(length(indic_weights) != 3 | !is.numeric(indic_weights)){
warning("Weights must be a numeric vector of length 3 for this model. Using the default (noise:1, population: 0.5, health: 0.5). Specify a weight of 0 if you want to ignore a dataset.", call. = FALSE)
indic_weights <- c(1, 0.5, 0.5)
}
)
# Are the weights specified correctly?
if (is.numeric(indic_weights)){
if (length(indic_weights) != 3){ # if only one weight specified we assign it to each layer
message("Weights should have 3 values for this model. Defaulting to weights of c(1, 0.5, 0.5).")
indic_weights <- c(1,1,1,1)
}
} else stop("Weights must be a numeric vector of length 3 for this model. Please refer to documentation. Specify a weight of 0 if you want to ignore a dataset.")
# Check the distance thresholds
if (!all(names(distances) == c("motorway", "dual", "primary", "railway", "airport")) | !is.numeric(distances)){
stop("Distance thresholds set incorrectly. Must be a named vector: make sure that names match documentation. Alternatively remove the argument to use default settings.")
}
# if mm is stored in list, combine all before proceeding
if (isTRUE(class(x) == "list")){
x <- do.call(rbind, x) %>% sf::st_as_sf()
x <- checkgeometry(x, "POLYGON")
# NOT using rbindlist here because only keeps the extent of the first tile
}
# Are any attributes missing?
if (any(!attributes %in% names(x))){
stop("Could not find attribute(s) in your basemap: ",
paste(attributes[which(!attributes %in% names(x))], collapse = ", "),
". Please add or rename the attributes.")
}
### Keep only fields we need
message("Preparing basemap and study area.")
x <- x[attributes]
# Create a buffer around the study area (to clip roads with some margin)
SAbuffer <- sf::st_buffer(sf::st_geometry(studyArea), 500) %>% sf::st_as_sf() %>% sf::st_transform(27700)
# Create raster with same properties as mastermap to use as a template
r <- raster::raster() # create empty raster
raster::crs(r) <- sp::CRS(SRS_string = "EPSG:27700") # hard-coding datum to preserve CRS
raster::extent(r) <- raster::extent(x) # set same extent as the basemap
raster::res(r) <- res # set resolution
################################# DISTANCE RASTER ##############################
if (indic_weights[1] > 0){ # only if user wants this indicator
# NOTE - Higher scores here mean an area is further from a noise source, so scores are inverted to mean that higher
# scores represent greater need for noise regulation
message("Importing noise sources (OS VectorMap) from EcoservR")
##### Importing OS VectorMap data #####
suppressWarnings({
airports <- readRDS(system.file("extdata/vectormap/airports.RDS", package = "ecoservR")) %>% sf::st_set_crs(27700)
railways <- readRDS(system.file("extdata/vectormap/railways.RDS", package = "ecoservR")) %>% sf::st_set_crs(27700)
roads <- readRDS(system.file("extdata/vectormap/roads.RDS", package = "ecoservR")) %>% sf::st_set_crs(27700)
})
#### Clip them to the buffered study area
airports <- sf::st_intersection(airports, SAbuffer)
railways <- sf::st_intersection(railways, SAbuffer)
roads <- sf::st_intersection(roads, SAbuffer)
##### Subsetting roads and buffering sources #####
# Create layers for noise sources
# (ONLY if there are occurrences in the study area, otherwise create empty raster)
message("Preparing VectorMap data for noise distance analysis")
motorways <- dplyr::filter(roads, classification == "Motorway") %>% # subset of roads that are motorways
sf::st_buffer(15) %>%
checkgeometry() # buffer and convert to polygon
dualcarriageways <- dplyr::filter(roads, classification == "A Road") %>% # subset of roads that are dual carriageways
sf::st_buffer(5) %>%
checkgeometry() # buffer and convert to polygon
roadsA <- dplyr::filter(roads, classification == "Minor Road") %>% # subset of A roads
sf::st_buffer(5) %>%
checkgeometry() # buffer and convert to polygon
airports <- airports %>%
sf::st_buffer(500) %>%
checkgeometry()
railways <- railways %>%
sf::st_buffer(5) %>%
checkgeometry()
##### 4 - Calculating noise source distances #####
message("Creating noise distance raster (indicator 1 of 3)")
# calculate distance to noise sources, writing resulting rasters to scratch folder:
# here the maxdist argument could become a reactive argument defined by thresholds from the defra noise dataset
dist_motor <- distance_from_source(motorways, r, "motorways", distances[["motorway"]])
dist_roadsA <- distance_from_source(roadsA, r, "roadsA", distances[["primary"]])
dist_dualcarriageways <- distance_from_source(dualcarriageways, r, "dualcarriageways", distances[["dual"]])
dist_rail <- distance_from_source(railways, r, "railways", distances[["railway"]])
dist_air <- distance_from_source(airports, r, "airports", distances[["airport"]])
### TO DO: REVISE DISTANCES SET ABOVE AND BELOW (THRESHOLDS FOR EACH SOURCE)
## This was actual distances but the score must be inversely proportional.
## From ecoservGIS:converted to log10 score and then to inverse score
# also added check because there might not be any features of a kind, and empty raster causes error in calculation
dist_motor <- if (raster::hasValues(dist_motor)){1 - (log10(dist_motor + 1) / log10(distances[["motorway"]]))} else {r}
dist_roadsA <- if (raster::hasValues(dist_roadsA)){1 - (log10(dist_roadsA + 1) / log10(distances[["primary"]]))} else {r}
dist_dualcarriageways <- if (raster::hasValues(dist_dualcarriageways)){1 - (log10(dist_dualcarriageways + 1) / log10(distances[["dual"]]))} else {r}
dist_rail <- if (raster::hasValues(dist_rail)){1 - (log10(dist_rail + 1) / log10(distances[["railway"]]))} else {r}
dist_air <- if (raster::hasValues(dist_air)){1 - (log10(dist_air + 1) / log10(distances[["airport"]]))} else {r}
##### 5 - Summing distance rasters #####
# Sum rasters:
# NB: have to use stack function (originally brick, caused a bug)
distances <- raster::stack(dist_motor, dist_roadsA, dist_dualcarriageways, dist_rail, dist_air) # stack rasters
rm(dist_motor, dist_roadsA, dist_dualcarriageways, dist_rail, dist_air)
if (terra::hasValues(distances)){
distancescore <- raster::calc(distances, fun = sum, na.rm = TRUE) #sum scores
# score could be slightly negative (precision), so forcing to zero
distancescore[distancescore < 0] <- 0
} else {
# if no noise sources at all in study area, set raster to 0
distancescore <- r
terra::values(distancescore) <- 0
}
# save raw indicator
if (indicators){
distancescore <- raster::writeRaster(distancescore,
filename = file.path(indicator_path,
paste(projectLog$title,runtitle, "noise_reg_dist_to_noise_indic.tif", sep="_")),
overwrite = TRUE)
message("Distance to noise sources indicator saved.")
}
rm(airports, motorways, railways, roads, roadsA, dualcarriageways)
} else {
distancescore <- r # if user wishes to omit this indicator we create empty raster
}
### POP AND/OR HEALTH
## If at least the pop and/or health indicators are selected (non-zero weights), create the houses
if (indic_weights[2] > 0 | indic_weights[3] > 0){
# Create the spatial points with socioeconomic data for each household
# We can't go straight from polygon to raster, because then a big house that spans two raster cells would
# have its population counted twice. By converting to points first, we make sure that each house is only represented once. If more than house house falls within a raster cell, their population gets summed.
houses <- dplyr::filter(x, HabCode_B == "J360") %>% # filter houses from basemap
checkgeometry("POLYGON") %>% # need to be single-part before converting to centro
sf::st_centroid() # get the centroid of each house
}
# Population raster -------------------------------------------------------
if (indic_weights[2] > 0){ # only if user wants this indicator
# Population size (sum) within local search distance
# NOTE - Higher population scores mean there is more societal need
message("INDICATOR 2 of 3: population density")
## Create the population raster by summing individual house populations in each cell
popscore <- raster::rasterize(houses, r, field = "housePop", fun = "sum", na.rm = TRUE,
filename = file.path(scratch, "airpurif_pop_temp"), overwrite = TRUE)
message("...calculating population in ", local, " m radius")
# Calculate the summarised local scores
popscore <- focalScore(popscore, radius = local, type = "sum")
# ## OLD FOCAL STATS
# # Create weight matrices based on the search radius for focal stats
# # (automatically considers the res of the raster to calculate distance)
# window_local <- raster::focalWeight(r, local, "circle") # search window for the longer dist
# window_local[window_local > 0] <- 1 # replacing weights by 1 (we want full sum, not mean)
#
#
# message("...calculating population in ", local, " m radius")
#
# # Calculate the summarised local scores
# popscore <- raster::focal(popscore, w = window_local, na.rm = TRUE, pad = TRUE)
# Save raw indicator if need be
if (indicators){
raster::writeRaster(popscore, file.path(indicator_path, # we define the path
paste(projectLog$title,runtitle, "noise_reg_pop_indic.tif", sep = "_")),
overwrite = TRUE)
}
} else {
popscore <- r
}
# Health indicator --------------------------------------------------------
if (indic_weights[3] > 0){ # only if user wants this indicator
# Health deprivation score (mean) within local search distance
# NOTE - Higher health scores mean an area is more deprived, i.e. there is more societal need for noise regulation
message("INDICATOR 3 of 3: health score")
# Create the health raster from the houses data we subsetted in previous indicator
healthscore <- raster::rasterize(houses, r, field = "health", fun = mean, na.rm = TRUE,
filename = file.path(scratch, "airpurif_health_temp"),
overwrite = TRUE)
rm(houses)
# Calculate the summarised local scores
healthscore <- focalScore(healthscore, radius = local, type = "mean")
### OLD FOCAL STATS
# healthscore <- raster::focal(healthscore, w = window_local,
# fun = function(x){mean(x[window_local != 0], na.rm=TRUE)
# # we select only the values in the circular window (discarding values where window is 0), and perform a mean calculation on this subset rather than the full rectangular extent of the matrix. This ensure the correct denominator while discarding NAs.
# },
# pad = TRUE) # pad edges with NAs so raster has full extent)
# Save raw indicator if need be
if (indicators){
raster::writeRaster(healthscore, file.path(indicator_path, # we define the path
paste(projectLog$title, runtitle, "noise_reg_health_indic.tif", sep = "_")),
overwrite = TRUE)
}
} else {
healthscore <- r
}
# RESCALE INDICATORS ------------------------------------------------------
# Each indicator gets transformed into a z-score (mean-centered and deviance-scaled), and then reprojected on a 1-100 scale
distancescore <- rescale_indicator(distancescore)
popscore <- rescale_indicator(popscore)
healthscore <- rescale_indicator(healthscore)
# COMBINE INDICATORS ------------------------------------------------------
# All indicators should now be between 0-1 and we can sum them.
message("Combining all indicators...")
final <- sum(
indic_weights[1]*distancescore,
indic_weights[2]*popscore,
indic_weights[3]*healthscore,
na.rm = TRUE # important otherwise only gives raster for where all 4 indicators have values
)
rm(distancescore, popscore, healthscore)
# Replace NA with 0
final[is.na(final)] <- 0
# Mask to study area
final <- raster::mask(final, studyArea)
message("Saving final scores.")
raster::writeRaster(final,
filename = file.path(save, paste(projectLog$title, runtitle,
"noise_reg_demand.tif",
sep="_")),
overwrite = TRUE)
## and a rescaled version
maxvalue <- max(raster::values(final), na.rm = TRUE)
raster::writeRaster(final/maxvalue *100,
filename = file.path(save, paste(projectLog$title, runtitle,
"noise_reg_demand_rescaled.tif", sep="_")),
overwrite = TRUE)
timeB <- Sys.time() # start time
# write performance to log
projectLog$performance[["dem_noise"]] <- as.numeric(difftime(
timeB, timeA, units="mins"
))
updateProjectLog(projectLog) # save revised log
# Delete all the stuff we don't need anymore
on.exit({
rm(r, final, maxvalue)
cleanUp(scratch)
message("Noise regulation demand model finished. Process took ", round(difftime(timeB, timeA, units = "mins"), digits = 1), " minutes. Please check output folder for your maps.")
})
return({
## returns the objects in the global environment
invisible({
projectLog <<- projectLog
})
})
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.