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

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(home2park)

devtools::load_all() # to knit manually w latest changes (not built/installed yet)
# vignette will be checked in R-CMD check & package build (should use library(home2park))

data(buildings_pop_sgp)
buildings_pop_sgp <- sf::st_transform(buildings_pop_sgp, sf::st_crs(32648)) # transform to projected crs

data(parks_sgp)
parks_sgp <- sf::st_transform(parks_sgp, sf::st_crs(32648)) # transform to projected crs

data(pop_sgp)
pop_sgp <- sf::st_transform(pop_sgp, sf::st_crs(32648)) # transform to projected crs

# convert buildings to points (centroids), then calculate distances to every park
m_dist <- buildings_pop_sgp %>%
  sf::st_centroid() %>%
  sf::st_distance(parks_sgp) %>% # euclidean distance
    units::set_units(NULL)

m_dist <- m_dist / 1000 # convert distances to km


# new column for the supply of park area
buildings_pop_sgp$area_supply <- recre_supply(park_attribute = parks_sgp$area, 
                                              dist_matrix = m_dist, 
                                              c = 0.302) # e.g. from Tu et al. (2020)

# supply to all residents per building
buildings_pop_sgp$area_supplytopop <- buildings_pop_sgp$area_supply * buildings_pop_sgp$popcount

rm(m_dist)

# calculate for year 2020
pop_2020 <- pop_sgp %>%
  dplyr::filter(year == 2020)

# append subzone info to buildings if building intersects with subzone  
buildings_subzones <- buildings_pop_sgp  %>%
  sf::st_join(pop_2020,
              join = sf::st_intersects,
              left = TRUE) %>%
  sf::st_set_geometry(NULL)


# summarise total/median/average supply value for all buildings within subzone
buildings_subzones <- buildings_subzones %>%
  dplyr::group_by(subzone_n) %>%
  dplyr::summarise(across(.cols = c(area_supply, area_supplytopop), 
                          .fns = sum, .names = "{.col}_sum"),
                   across(.cols = c(area_supply, area_supplytopop), 
                          .fns = median, .names = "{.col}_median"),
                   across(.cols = c(area_supply, area_supplytopop), 
                          .fns = mean, .names = "{.col}_mean")) 


# join information to pop_2020, calculate supply per capita
pop_2020 <- pop_2020 %>%
  dplyr::left_join(buildings_subzones) %>%
  dplyr::mutate(area_supplyperpop = area_supplytopop_sum / pop_count * 1e-6) %>% # convert to km2
  dplyr::mutate(area_supplyperpop = ifelse(is.infinite(area_supplyperpop), NA, area_supplyperpop)) # make infinite NA

sf::st_agr(pop_2020) = "constant"
sf::st_agr(parks_sgp) = "constant"
# https://github.com/r-spatial/sf/issues/406

subzones_parks <- sf::st_intersection(pop_2020, parks_sgp) # subset census unit polygons to those that intersect parks, append park info
subzones_parks$parkarea_m2 <- sf::st_area(subzones_parks) # calc area of each polygon


# calculate total park area per census unit
subzones_parks <- subzones_parks %>%
  dplyr::group_by(subzone_n) %>%
  dplyr::summarise(parkarea_m2 = as.numeric(sum(parkarea_m2))) %>%
  sf::st_drop_geometry()


# join information to pop_2020, calculate park area per capita
pop_2020 <- pop_2020 %>%
  dplyr::left_join(subzones_parks) %>%
  dplyr::mutate(parkarea_m2 = ifelse(is.na(parkarea_m2), 0, parkarea_m2)) %>% # make NA 0
  dplyr::mutate(parkperpop_m2 = parkarea_m2 / pop_count * 1e-6) %>% # convert to km2
  dplyr::mutate(parkperpop_m2 = ifelse(is.infinite(parkperpop_m2), # make infinite NA
                                NA, parkperpop_m2)) %>%
  dplyr::mutate(parkperpop_m2 = ifelse(is.nan(parkperpop_m2), # make NaN NA
                                NA, parkperpop_m2))

rm(buildings_subzones, subzones_parks)

The following interactive map provides a visual comparison between metrics for park area provision calculated in the 'Get started' vignette. Toggle the map layers to view the per capita provision of park area summarised per region (subzone).

# convert buildings to points (centroids) for plotting
buildings_pop_sgp <- buildings_pop_sgp %>% 
  sf::st_centroid() %>%
  dplyr::mutate(across(everything(), as.vector)) # remove attributes from columns for plotting


# for area related vars, convert m2 to km2
buildings_pop_sgp <- buildings_pop_sgp %>% 
  dplyr::mutate(across(.cols = contains("area"), 
              .fns = function(x) x*1e-6))


# random sampling of buildings
# set.seed(123) 
# # some processing to get proper range for color scale later
# # get max values per building (for particular round)
# buildings_max <- buildings_pop_sgp %>%
#   tidyr::pivot_longer(cols = c("area_supply", "area_supplytopop"),
#                       names_to = "supply", 
#                       values_to = "value") %>%
#      dplyr::group_by(supply) %>%
#      dplyr::slice(which.max(value)) %>%
#   tidyr::pivot_wider(names_from = "supply", values_from = "value")
# 
# # subset random sample (already done with example dataset)
# buildings_pop_sgp <-
#   dplyr::slice_sample(buildings_pop_sgp,
#                       n = nrow(buildings_pop_sgp)/20)
# 
# # append max buildings to the random subset 
# buildings_pop_sgp <- buildings_pop_sgp %>%
#   dplyr::bind_rows(buildings_max)





tmap::tmap_mode("view")
tmap::tmap_options(check.and.fix = TRUE)

tm <- 
  tmap::tm_basemap(c("Esri.WorldGrayCanvas", "CartoDB.DarkMatter", "OpenStreetMap")) +

  # parks  
  tmap::tm_shape(parks_sgp %>% dplyr::select(id, name, area)) +
    tmap::tm_polygons(group = "Parks",
                      col = "#33a02c",
                      alpha = 0.6,
                      border.col = "grey50",  border.alpha = 0.5) +

  # buildings
  tmap::tm_shape(buildings_pop_sgp) +
    tmap::tm_dots(title = "Buildings: Supply of park area (km<sup>2</sup>)",
                  group = "Buildings: Supply of park area",
                  col = "area_supply", border.col = "transparent",
                  palette = viridis::viridis(5),
                  style = "quantile",
                  size = 0.01,
                  alpha = 0.7,
                  interactive = FALSE, # some glitch in values of hover text/popup
                  popup.vars = NULL,
                  showNA = FALSE) +

  tmap::tm_shape(buildings_pop_sgp) +
    tmap::tm_dots(title = "All building residents: Supply of park area (km<sup>2</sup>)",
                  group = "Buildings: Supply of park area to residents",
                  col = "area_supplytopop", border.col = "transparent",
                  palette = viridis::viridis(5),
                  style = "quantile",
                  size = 0.01,
                  alpha = 0.7,
                  interactive = FALSE, # some glitch in values of hover text/popup
                  popup.vars = NULL,
                  showNA = FALSE) +

  # census units
  tmap::tm_shape(pop_2020) +
    tmap::tm_polygons(title = "Supply of park area per capita (km<sup>2</sup>)",
                      group = "Subzones: Per capita supply of park area",
                      col = "area_supplyperpop", 
                      palette = "Greens", alpha = 0.7,
                      style = "quantile",
                      border.col = "white", border.alpha = 0.5, lwd = 1) +
  tmap::tm_shape(pop_2020) +
    tmap::tm_polygons(title = "Park area per capita (km<sup>2</sup>)",
                      group = "Subzones: Per capita park area (conventional)",
                      col = "parkperpop_m2",
                      palette = "Greens", alpha = 0.7,
                      style = "quantile",
                      border.col = "white", border.alpha = 0.5, lwd = 1)


# Pipe the tmap object into tmap_leaflet() to create a leaflet widget,
# so that we can use leaflet::hideGroup().
tm %>% 
  tmap::tmap_leaflet() %>%
  leaflet::hideGroup("Buildings: Supply of park area") %>%
  leaflet::hideGroup("Subzones: Per capita supply of park area") %>% 
  leaflet::hideGroup("Subzones: Per capita park area (conventional)")

rm(tm, buildings_pop_sgp, parks_sgp, pop_sgp, pop_2020)

Data sources

Singapore census data from the Department of Statistics Singapore. Released under the terms of the Singapore Open Data Licence version 1.0.
Singapore subzone polygons from the Singapore Master Plan Subzones. Released under the terms of the Singapore Open Data Licence version 1.0.
Singapore Master Plan Land Use Zones for the years 2014 and 2019. Released under the terms of the Singapore Open Data License.
Building polygons derived from map data copyrighted OpenStreetMap contributors and available from https://www.openstreetmap.org. Released under the terms of the ODbL License.
Park polygons and summarised attributes (trails, playgrounds) derived from map data copyrighted OpenStreetMap contributors and available from https://www.openstreetmap.org. Released under the terms of the ODbL License.