Visualizing Highways and Healthcare Facilities

Analyzing the integration of healthcare facilities and transportation networks - using power of {osmdata}, {sf} and {ggspatial}

Data Visualization

Maps

Public Health

{gt}

Gecomputation

Author

Aditya Dahiya

Published

November 24, 2024

Introduction

This project visualizes the distribution of health facilities in Haryana, India, alongside the national highways, using data from OpenStreetMap. The primary aim is to map health facilities categorized as small, medium, or large and assess their proximity to highways (within 200 meters or farther).

Process:

Data Collection: Using the osmdata package, bounding boxes for Haryana, Delhi and later India were defined. Data on highways and health facilities were downloaded, pre-processed, and saved for subsequent analysis.
Spatial Processing: The sf package (Pebesma and Bivand 2023) facilitated spatial operations, such as clipping datasets to Delhi’s, Haryana’s and India’s boundaries, transforming geometry types (e.g., polygons to points on the map using st_centroid()), and, Most importantly, calculating distances between health facilities and highways using st_is_within_distance()with techniques from the book Geocomputation with R (Chapter 4).
Visualization: A detailed map was created with ggplot2, incorporating highways, health facilities, and annotations for scale and direction using ggspatial. Facilities were color-coded (near or far from highways) and differentiated by size and shape.
Inset Charts: Additional insights were provided through inset charts showing the number of health facilities by district in Haryana, and the percentage of facilities near highways, categorized by size, across India.

Key Insights:

Proximity Analysis: The project quantifies the percentage of health facilities near highways, highlighting disparities in accessibility.
Facility Distribution: The data shows varying distributions of small, medium, and large facilities across districts, providing insights into healthcare accessibility.
Comprehensive Visualization: The use of patchwork to combine plots and gt for tabular data enhances the project’s visual and analytical depth.

The final map, saved as an image using ggsave, effectively communicates the spatial relationship between health infrastructure and transportation networks. This visualization can inform infrastructure planning and policy decisions.

Loading Libraries

# Data Import and Wrangling Tools
library(tidyverse)            # All things tidy
library(here)                 # File locations and paths

# Final plot tools
library(scales)               # Nice Scales for ggplot2
library(fontawesome)          # Icons display in ggplot2
library(ggtext)               # Markdown text support for ggplot2
library(showtext)             # Display fonts in ggplot2
library(colorspace)           # Lighten and Darken colours
library(patchwork)            # Combining plots
library(gt)                   # Beautiful Tables in R

# Mapping tools
library(sf)                   # All spatial objects in R
library(osmdata)              # Getting Open Street Maps data

Part 1: Getting Raw Data

This code demonstrates how to query and store OpenStreetMap (OSM) data for health amenities and highways in Haryana and Delhi using the osmdata package ( et al. 2017). Initially, bounding boxes for Haryana (hy_bbox) and Delhi (del_bbox) are obtained using osmdata::getbb(). For each region, two categories of health facilities are retrieved: smaller establishments like clinics, dentists, and pharmacies (health_small) and larger institutions like hospitals and nursing homes (health_large). These datasets are fetched using the opq() function to create an Overpass API query, combined with add_osm_feature() to specify the type of amenities, and are then converted into simple feature (sf) objects via osmdata_sf().

Similarly, highways are categorized into major roads (motorways and trunks) and secondary roads (primary and secondary), and their data is retrieved and converted into sf objects. To ensure efficient future use, all datasets are temporarily saved as RDS files using saveRDS(). This structured approach ensures the data is readily accessible for visualization or analysis. For more on querying OSM data, refer to the osmdata package documentation.

Code

hy_bbox <- osmdata::getbb("Haryana")

del_bbox <- osmdata::getbb("Delhi")

######## Haryana Health
health_small <- opq(bbox = hy_bbox) |>
  add_osm_feature(
    key = "amenity",
    value = c("clinic", "dentist", "doctors", "pharmacy")
  ) |>
  osmdata_sf()

object.size(health_small) |> print(units = "Mb")

health_large <- opq(bbox = hy_bbox) |>
  add_osm_feature(
    key = "amenity",
    value = c("hospital", "nursing_home", "social_facility")
  ) |>
  osmdata_sf()

object.size(health_large) |> print(units = "Mb")

# Temporary Saving the Data for easy use in travel
saveRDS(
  object = health_large,
  file = "haryana_health_large.rds"
)

saveRDS(
  object = health_small,
  file = "haryana_health_small.rds"
)

# Haryana Highways ---------------------
roads_1 <- opq(bbox = hy_bbox) |> 
  add_osm_feature(
    key = "highway", 
    value = c("motorway", "trunk", "motorway_link", "trunk_link")
  ) |> 
  osmdata_sf()

object.size(roads_1) |> print(units = "Mb")

roads_2 <- opq(bbox = hy_bbox) |> 
  add_osm_feature(
    key = "highway", 
    value = c("primary", "secondary")
  ) |> 
  osmdata_sf()

object.size(roads_2) |> print(units = "Mb")

# Temporary Saving the Data for easy use in travel
saveRDS(
  object = roads_1,
  file = "haryana_roads_1.rds"
)

saveRDS(
  object = roads_2,
  file = "haryana_roads_2.rds"
)

############ Delhi Health
health_small <- opq(bbox = del_bbox) |>
  add_osm_feature(
    key = "amenity",
    value = c("clinic", "dentist", "doctors", "pharmacy")
  ) |>
  osmdata_sf()

object.size(health_small) |> print(units = "Mb")

health_large <- opq(bbox = del_bbox) |>
  add_osm_feature(
    key = "amenity",
    value = c("hospital", "nursing_home", "social_facility")
  ) |>
  osmdata_sf()

object.size(health_large) |> print(units = "Mb")

# Temporary Saving the Data for easy use in travel
saveRDS(
  object = health_large,
  file = "delhi_health_large.rds"
)

saveRDS(
  object = health_small,
  file = "delhi_health_small.rds"
)

# Delhi Highways ---------------------
roads_1 <- opq(bbox = del_bbox) |> 
  add_osm_feature(
    key = "highway", 
    value = c("motorway", "trunk", "motorway_link", "trunk_link")
  ) |> 
  osmdata_sf()

object.size(roads_1) |> print(units = "Mb")

roads_2 <- opq(bbox = del_bbox) |> 
  add_osm_feature(
    key = "highway", 
    value = c("primary", "secondary")
  ) |> 
  osmdata_sf()

object.size(roads_2) |> print(units = "Mb")

# Temporary Saving the Data for easy use in travel
saveRDS(
  object = roads_1,
  file = "delhi_roads_1.rds"
)

saveRDS(
  object = roads_2,
  file = "delhi_roads_2.rds"
)

Part 2: Get massive data for India

This code focuses on querying and saving OpenStreetMap (OSM) data for health amenities and highways across entire India. To ensure efficient handling of such a large geographical scope, the data retrieval process is divided into smaller datasets. Health facilities are then split into two categories: smaller facilities like clinics, dentists, and pharmacies (health_small) and larger institutions such as hospitals, nursing homes, and social facilities (health_large). For highways, the script retrieves major road types, including motorways and trunk roads (roads_1) only. While a placeholder for secondary roads exists (roads_2), it is commented out to optimize download times or avoid potential interruptions. To facilitate future use, all retrieved datasets are saved as RDS files using saveRDS().

Code

###########################################################
###############           India           #################
###########################################################

in_bbox <- osmdata::getbb("India")

######## India Health -----------------------------------------

# Note: I have split downloading health facilties into larger and
#       Smaller ones to save on downloading times, and prevent 
#       breakage during downloading such large datasets

health_small <- opq(bbox = in_bbox) |>
  add_osm_feature(
    key = "amenity",
    value = c("clinic", "dentist", "doctors", "pharmacy")
  ) |>
  osmdata_sf()

object.size(health_small) |> print(units = "Mb")

# Temporarily Saving the Data for easy use in travel / reload
saveRDS(
  object = health_small,
  file = "india_health_small.rds"
)

health_large <- opq(bbox = in_bbox) |>
  add_osm_feature(
    key = "amenity",
    value = c("hospital", "nursing_home")
  ) |>
  osmdata_sf()

object.size(health_large) |> print(units = "Mb")

# Temporarily Saving the Data for easy use in travel / reload
saveRDS(
  object = health_large,
  file = "india_health_large.rds"
)

########### India Highways -----------------------------------
# Troubleshoot: Large data, hence need to increase time-out
# Credits: https://github.com/ropensci/osmdata/issues/200
roads_1 <- opq(bbox = in_bbox, timeout = 500) |> 
  add_osm_feature(
    key = "highway", 
    value = c("motorway")
  ) |> 
  osmdata_sf()

object.size(roads_1) |> print(units = "Mb")

saveRDS(
  object = roads_1,
  file = "india_roads_1.rds"
)

roads_1_5 <- opq(bbox = in_bbox, timeout = 100) |> 
  add_osm_feature(
    key = "highway", 
    value = c("trunk")
  ) |> 
  osmdata_sf()

object.size(roads_1) |> print(units = "Mb")

saveRDS(
  object = roads_1_5,
  file = "india_roads_1_5.rds"
)
# 
# roads_2 <- opq(bbox = in_bbox) |> 
#   add_osm_feature(
#     key = "highway", 
#     value = c("primary", "secondary")
#   ) |> 
#   osmdata_sf()
# 
# object.size(roads_2) |> print(units = "Mb")
# 
# # Temporary Saving the Data for easy use in travel
# 
# saveRDS(
#   object = roads_2,
#   file = "india_roads_2.rds"
# )

Part 3: Trial Code with Delhi’s smaller dataset

The process begins by loading the saved datasets for health facilities (health_small and health_large) and highways (roads_1 and roads_2) using readRDS(). Each dataset is further refined to include only the relevant geographic data (osm_lines for highways and points or centroids for health facilities).

A shapefile of Delhi’s boundary is imported and cleaned, ensuring that the spatial coordinate reference systems match for all datasets. Highways and health facilities are filtered to include only those within Delhi’s boundary using st_intersection().

For health facilities, geometries (e.g., polygons) are converted to centroids, ensuring consistency in spatial representation. Facilities are classified based on size (e.g., “small”, “medium”, or “large”) and evaluated for their proximity to highways using st_is_within_distance(). A distance threshold of 200 meters is applied to categorize facilities as "near_highway" or "far_from_highway".

The results are summarized in tabular format for both small and large health facilities, showing the distribution of facilities based on proximity to highways. Finally, the data is visualized in Figure 1 using ggplot2, with:

Highways (roads_1 and roads_2) displayed as lines of varying weights and colors.
Health facilities plotted as points, styled by their proximity to highways.

Code

# Get pre-saved data (Trial Phase)
health_small <- readRDS(file = "delhi_health_small.rds")
health_large <- readRDS(file = "delhi_health_large.rds")
roads_1 <- readRDS(file = "delhi_roads_1.rds")
roads_2 <- readRDS(file = "delhi_roads_2.rds")

# Get only the relevant sf obejct out of each
roads_1 <- roads_1$osm_lines
roads_2 <- roads_2$osm_lines

# Get a Map of the State / City and keep only the highways and facilties within that State / City
region_map <- sf::read_sf(
  here::here("data", "india_map", "India_State_Boundary.shp")
) |> 
  janitor::clean_names() |> 
  filter(state_name == "Delhi") |> 
  st_transform(crs = st_crs(roads_1))

# ggplot() +
#   geom_sf(data = region_map) +
#   geom_sf(data = roads_1)
# 
# ggplot() +
#   geom_sf(data = region_map) +
#   geom_sf(data = roads_1 |> st_intersection(region_map))

# Select Roads and Facilties that are within the selected region
roads_1 <- roads_1 |> st_intersection(region_map)
roads_2 <- roads_2 |> st_intersection(region_map)

# Health facilities are as points, polygons and multi-polygons. 
# Lets convert them all into points, since on a large map 
# like a city or state, all will anyways appear as a point

# Exploration of data available in the health facilities
# health_large$osm_points |> 
#   st_drop_geometry() |> 
#   select(name, beds, building, `building:levels`,highway, short_name) |> 
#   visdat::vis_miss()
#   names()

# Do for health_large -------------------------------------------
temp1 <- health_large$osm_points |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "medium")

temp2 <- health_large$osm_polygons |> 
  st_centroid() |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "large")

temp3 <- health_large$osm_multipolygons |> 
  st_centroid() |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "large")

health_large <- bind_rows(
  temp1,
  temp2,
  temp3
)

rm(temp1, temp2, temp3)

health_large <- health_large |> 
  st_intersection(region_map)

health_large <- health_large |> 
  mutate(
    near_highway = health_large |> 
      
      # Checking which points are within a specific distance of highway
      st_is_within_distance(
        y = roads_1,
        dist = 200,
        sparse = F
      ) |> 
      as_tibble() |> 
      mutate(
        any_true = rowSums(across(everything(), as.logical)) > 0
      ) |> 
      pull(any_true)
  ) |> 
  mutate(
    near_highway = if_else(
      near_highway,
      "near_highway",
      "far_from_highway"
    )
  )

# A Table of near and away from highway facilties for 
# Large and Medium Health Facilties
table_large <- health_large |> 
  st_drop_geometry() |> 
  count(facility_type, near_highway) |> 
  pivot_wider(
    id_cols = facility_type,
    names_from = near_highway,
    values_from = n
  ) |> 
  mutate(
    perc_on_highway = near_highway / (near_highway + far_from_highway)
  )

# Do for health_small --------------------------------------------
temp1 <- health_small$osm_points |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "small")

temp2 <- health_small$osm_polygons |> 
  st_centroid() |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "medium")

# temp3 <- health_small$osm_multipolygons |> 
#   st_centroid() |> 
#   select(name, geometry) |> 
#   mutate(facility_type = "medium")

health_small <- bind_rows(
  temp1,
  temp2
)

rm(temp1, temp2)

health_small <- health_small |> 
  st_intersection(region_map)

health_small <- health_small |> 
  mutate(
    near_highway = health_small |> 
      
      # Checking which points are within a specific distance of highway
      st_is_within_distance(
        y = roads_1,
        dist = 200,
        sparse = F
      ) |> 
      as_tibble() |> 
      mutate(
        any_true = rowSums(across(everything(), as.logical)) > 0
      ) |> 
      pull(any_true)
  ) |> 
  mutate(
    near_highway = if_else(
      near_highway,
      "near_highway",
      "far_from_highway"
    )
  )


# A Table of near and away from highway facilties for 
# Large and Medium Health Facilties
table_small <- health_small |> 
  st_drop_geometry() |> 
  count(facility_type, near_highway) |> 
  pivot_wider(
    id_cols = facility_type,
    names_from = near_highway,
    values_from = n
  ) |> 
  mutate(
    perc_on_highway = near_highway / (near_highway + far_from_highway)
  )

# The final table product -------------------------------------------
show_table <- bind_rows(table_large, table_small) |> 
  group_by(facility_type) |> 
  summarise(
    far_from_highway = sum(far_from_highway),
    near_highway = sum(near_highway)
  ) |> 
  mutate(
    perc_on_highway = near_highway / (near_highway + far_from_highway)
  )

# The plot

g <- ggplot() +
  geom_sf(
    data = region_map,
    linewidth  = 0.5,
    colour = "grey20",
    fill = "transparent"
  ) +
  geom_sf(
    data = roads_1, 
    linewidth = 0.5, 
    alpha = 0.8,
    colour = "darkgrey"
  ) +
  geom_sf(
    data = roads_2, 
    linewidth = 0.3, 
    alpha = 0.4,
    colour = "lightgrey"
  ) +
  geom_sf(
    data = health_large,
    mapping = aes(
      colour = near_highway,
      size = near_highway
    ),
    alpha = 0.8
  ) +
  scale_size_manual(
    values = c(0.5, 1.5)
  ) +
  labs(
    title = "Delhi: Highways and Health Facilties"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom"
  )
  
ggsave(
  plot = g,
  filename = here("projects", "images", "highways_heath_1.png"),
  height = 1500,
  width = 1500,
  units = "px"
)

Part 4: Similar analysis for Haryana state

This R code analyzes the proximity of health facilities to highways in Haryana. It loads pre-saved datasets for small and large health facilities and roads, filters them by Haryana’s boundary, and converts all facilities into point geometries for simplicity. The proximity of facilities to highways is assessed using st_is_within_distance from the {sf} package. The code then generates summary tables categorizing facilities by type and proximity to highways. A map is created using {ggplot2} and {ggspatial}, showing highways and health facilities with annotations for scale and orientation. The final map is shown in Figure 2.

Code

# Get pre-saved data (Trial Phase - use pre-saved data to 
# save laoding times)
# Note: I have split downloading health facilties into larger and
#       Smaller ones to save on downloading times, and prevent 
#       breakage during downloading such large datasets
health_small <- readRDS(file = "haryana_health_small.rds")
health_large <- readRDS(file = "haryana_health_large.rds")
roads_1 <- readRDS(file = "haryana_roads_1.rds")
roads_2 <- readRDS(file = "haryana_roads_2.rds")

# Get only the relevant sf object out of each downlaoded data-set
# For roads, only "LINES" are relevant sf feature.
roads_1 <- roads_1$osm_lines
roads_2 <- roads_2$osm_lines

# Get a Map of the State / City and keep only the highways and 
# facilties within that State / City for plotting (instead of entire)
# bounding box
region_map <- sf::read_sf(
  here::here("data", "india_map", "India_State_Boundary.shp")
) |> 
  janitor::clean_names() |> 
  filter(state_name == "Haryana") |> 
  st_transform(crs = st_crs(roads_1))

# Checking the results
# ggplot() +
#   geom_sf(data = region_map) +
#   geom_sf(data = roads_1)
# 
# ggplot() +
#   geom_sf(data = region_map) +
#   geom_sf(data = roads_1 |> st_intersection(region_map))

# Select Roads and Facilties that are within the selected region
roads_1 <- roads_1 |> st_intersection(region_map)
roads_2 <- roads_2 |> st_intersection(region_map)

# Health facilities are as points, polygons and multi-polygons. 
# Lets convert them all into points, since on a large map 
# like a city or state, all will anyways appear as a point

# -------------------------------------------------------------------
# Extracting the Larger Health Facilities in a nicely formatted clean
# sf object, to make plotting easier subsequently--------------------
temp1 <- health_large$osm_points |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "medium")

temp2 <- health_large$osm_polygons |> 
  st_centroid() |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "large")

temp3 <- health_large$osm_multipolygons |> 
  st_centroid() |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "large")

health_large <- bind_rows(
  temp1,
  temp2,
  temp3
)
# Clean the workspace
rm(temp1, temp2, temp3)

# Retain only the health facilities that are within our region
health_large <- health_large |> 
  st_intersection(region_map)

# Now the MAGIC of {sf}: mark the facilities that are near highways
health_large <- health_large |> 
  mutate(
    near_highway = health_large |> 
      
      # Checking which points are within a specific distance of highway
      # The awesome power & speed of sf::st_is_within_distance() 
      st_is_within_distance(
        y = roads_1,
        dist = 200,
        sparse = F
      ) |> 
      as_tibble() |> 
      # Keep value TRUE if facility is near ANY highway
      mutate(
        any_true = rowSums(across(everything(), as.logical)) > 0
      ) |> 
      pull(any_true)
  ) |> 
  mutate(
    near_highway = if_else(
      near_highway,
      "near_highway",
      "far_from_highway"
    )
  )

# An Output Table of near and away from highway facilties for 
# Large and Medium Health Facilties
table_large <- health_large |> 
  st_drop_geometry() |> 
  count(facility_type, near_highway) |> 
  pivot_wider(
    id_cols = facility_type,
    names_from = near_highway,
    values_from = n
  ) |> 
  mutate(
    perc_on_highway = near_highway / (near_highway + far_from_highway)
  )

# Do the same process for snmaller health facilities ----------------
temp1 <- health_small$osm_points |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "small")

temp2 <- health_small$osm_polygons |> 
  st_centroid() |> 
  select(name, geometry) |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "medium")

health_small <- bind_rows(
  temp1,
  temp2
)

rm(temp1, temp2)

health_small <- health_small |> 
  st_intersection(region_map)

health_small <- health_small |> 
  mutate(
    near_highway = health_small |> 
      
      # Checking which points are within a specific distance of highway
      st_is_within_distance(
        y = roads_1,
        dist = 200,
        sparse = F
      ) |> 
      as_tibble() |> 
      mutate(
        any_true = rowSums(across(everything(), as.logical)) > 0
      ) |> 
      pull(any_true)
  ) |> 
  mutate(
    near_highway = if_else(
      near_highway,
      "near_highway",
      "far_from_highway"
    )
  )

# An output Table of near and away from highway facilties for 
# Large and Medium Health Facilties
table_small <- health_small |> 
  st_drop_geometry() |> 
  count(facility_type, near_highway) |> 
  pivot_wider(
    id_cols = facility_type,
    names_from = near_highway,
    values_from = n
  ) |> 
  mutate(
    perc_on_highway = near_highway / (near_highway + far_from_highway)
  )


# The final table product -------------------------------------------
show_table <- bind_rows(table_large, table_small) |> 
  group_by(facility_type) |> 
  summarise(
    far_from_highway = sum(far_from_highway),
    near_highway = sum(near_highway)
  ) |> 
  mutate(
    perc_on_highway = near_highway / (near_highway + far_from_highway)
  ) |> 
  mutate(
    facility_type = case_when(
      facility_type == "large" ~ "Large establishments: Hospitals, Nursing Homes, Institutions",
      facility_type == "medium" ~ "Medium: Clinics, Doctors's Practices",
      facility_type == "small" ~ "Small facilties: Dentists, Pharmacies",
      .default = NA
    )
  )

# Combining health_large and health_small to map shape to facility size
health_combined <- bind_rows(
  health_large,
  health_small
  ) |> 
  mutate(
    facility_type = fct(
      facility_type,
      levels = c("small", "medium", "large")
    )
  )

rm(health_large, health_small)

# The plot
g <- ggplot() +
  geom_sf(
    data = region_map,
    linewidth  = 0.5,
    colour = "grey20",
    fill = "transparent"
  ) +
  geom_sf(
    data = roads_1, 
    linewidth = 0.5, 
    alpha = 0.8,
    colour = "darkgrey"
  ) +
  geom_sf(
    data = roads_2,
    linewidth = 0.3,
    alpha = 0.4,
    colour = "lightgrey"
  ) +
  geom_sf(
    data = health_combined,
    mapping = aes(
      colour = near_highway,
      size = facility_type,
      shape = facility_type
    ),
    alpha = 0.7
  ) +
  scale_size_manual(
    values = c(0.2, 1, 2)
  ) +
  scale_shape_manual(
    values = c(20, 2, 0)
  ) +
  scale_colour_manual(
    values = c("lightblue", "red")
  ) +
  # Add North-Arrow and Scale
  ggspatial::annotation_north_arrow(
    style = ggspatial::north_arrow_orienteering(
      line_col = "grey20",
      fill = c("grey90", "grey10")
    ),
    location = "tl"
  ) +
  ggspatial::annotation_scale(
    bar_cols = c("grey10", "white"),
    location = "bl"
  ) +
  labs(
    title = "Haryana: Highways and Health Facilties",
    shape = NULL,
    size = NULL,
    colour = NULL
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "vertical"
  )

ggsave(
  plot = g,
  filename = here("projects", "images", "highways_heath_2.png"),
  height = 2500,
  width = 2000,
  units = "px"
)

Part 5: Try code for a {gt} table and an inset plot

This R code generates a table and a plot to analyze the proximity of health facilities to highways in Haryana. First, a summary table is created using {gt} with bar plots showing the percentage of health facilities located near highways. The table is formatted using gtExtras::gt_theme_538() for styling. Then, a plot is created with {ggplot2}, visualizing the percentage of health facilities located near highways by facility type (small, medium, large) and state. The data is transformed using {tidyverse} functions such as pivot_wider and count, and the plot is faceted by facility type for clarity.

Code

insert_gt_table <- show_table |> 
  mutate(facility_type = str_to_title(facility_type)) |> 
  rename(percentage_on_highway = perc_on_highway) |> 
  gt() |> 
  cols_label_with(fn = snakecase::to_title_case) |> 
  # cols_width(
  #   facility_type ~ px(100),
  #   far_from_highway ~ px(20),
  #   near_highway ~ px(20),
  #   percentage_on_highway ~ px(50)
  # ) |> 
  gtExtras::gt_plt_bar(
    percentage_on_highway,
    labels = TRUE,
    color = "darkgreen",
    scale_type = "percent"
  ) |> 
  gtExtras::gt_theme_538()

plot_data_df <- health_combined |> 
  st_drop_geometry() |> 
  as_tibble() |> 
  count(state_name, facility_type, near_highway) |> 
  pivot_wider(
    id_cols = c(state_name, facility_type),
    names_from = near_highway,
    values_from = n
  ) |> 
  mutate(
    total = far_from_highway + near_highway,
    perc_on_hwy = near_highway / total,
    .keep = "unused"
  )

strip_labels <- c(
      "Small facilties: Dentists, Pharmacies",
      "Medium: Clinics, Doctors's Practices",
      "Large: Hospitals, Nursing Homes, Institutions"
    ) |> 
  str_wrap(30)
names(strip_labels) <- c("small", "medium", "large")

inset_plot <- plot_data_df |> 
  ggplot(
    mapping = aes(
      y = state_name,
      x = perc_on_hwy
    )
  ) +
  geom_col(
    alpha = 0.5
  ) +
  geom_text(
    mapping = aes(
      label = paste0("(", number(total, 
                                 scale = 1, 
                                 accuracy = 1,
                                 big.mark = ","),
                     ")"),
      x = 0
    ),
    nudge_x = 0.01,
    hjust = 0
  ) +
  geom_text(
    mapping = aes(
      label = paste0(round(perc_on_hwy * 100, 2), "%")
    ),
    nudge_x = 0.005,
    hjust = 0
  ) +
  labs(
    y = NULL,
    x = "Health facilities located on the Highway"
  ) +
  facet_wrap(
    ~ facility_type,
    labeller = labeller(
      facility_type = strip_labels
    )
  ) +
  scale_x_continuous(
    expand = expansion(c(0, 0.25)),
    labels = scales::label_percent()
  ) +
  theme_minimal(
    base_size = 10
  ) +
  theme()

Part 6: A final product map for Haryana

The code creates a detailed Figure 3 for Haryana, showcasing highways and healthcare facilities. The map displays highways in different colours and health facilities marked by size and shape based on their type (small, medium, large) and proximity to highways (within 200 meters). The final map includes district boundaries, highways, and health facilities, with color indicating proximity to highways. A secondary plot summarizes healthcare facilities by district and type.

Code

# Basics of Final Visualization: Fonts
font_add_google("Rouge Script", "title_font")
font_add_google("El Messiri", "body_font")
font_add_google("Saira Extra Condensed", "caption_font")
showtext_auto()

text_col <- "grey10"
text_hil <- "grey20"
bts <- 150
bg_col <- "white"

# Caption stuff for the plot
sysfonts::font_add(
  family = "Font Awesome 6 Brands",
  regular = here::here("docs", "Font Awesome 6 Brands-Regular-400.otf")
)
github <- "&#xf09b"
github_username <- "aditya-dahiya"
xtwitter <- "&#xe61b"
xtwitter_username <- "@adityadahiyaias"
social_caption_1 <- glue::glue("<span style='font-family:\"Font Awesome 6 Brands\";'>{github};</span> <span style='color:{text_hil}'>{github_username}  </span>")
social_caption_2 <- glue::glue("<span style='font-family:\"Font Awesome 6 Brands\";'>{xtwitter};</span> <span style='color:{text_hil}'>{xtwitter_username}</span>")

plot_title <- "Haryana: Highways and Health-Care"

plot_subtitle <- "The purple lines depict highways in Haryana. Each dot represents a Health-Care Facility, with shape & size  depicting facility's level, and colour depicting proximity to a National Highway (within 200 metres)"

plot_caption <- paste0(
  "**Data:** Open Street Maps; Survey of India", 
  " |  **Code:** ", 
  social_caption_1, 
  " |  **Graphics:** ", 
  social_caption_2
  )

rm(github, github_username, xtwitter, 
   xtwitter_username, social_caption_1, 
   social_caption_2)

# Get pre-saved data (Trial Phase - use pre-saved data to 
# save laoding times)
# Note: I have split downloading health facilties into larger and
#       Smaller ones to save on downloading times, and prevent 
#       breakage during downloading such large datasets
health_small <- readRDS(file = "haryana_health_small.rds")
health_large <- readRDS(file = "haryana_health_large.rds")
roads_1 <- readRDS(file = "haryana_roads_1.rds")
roads_2 <- readRDS(file = "haryana_roads_2.rds")

# Get only the relevant sf object out of each downlaoded data-set
# For roads, only "LINES" are relevant sf feature.
roads_1 <- roads_1$osm_lines
roads_2 <- roads_2$osm_lines

# Get a Map of the State / City and keep only the highways and 
# facilties within that State / City for plotting (instead of entire)
# bounding box
region_map <- sf::read_sf(
  here::here("data", "india_map", "India_State_Boundary.shp")
) |> 
  janitor::clean_names() |> 
  filter(state_name == "Haryana") |> 
  st_transform(crs = st_crs(roads_1))


haryana_districts_map <- sf::read_sf(
  here::here("data", "haryana_map", "HARYANA_DISTRICT_BDY.shp")
) |> 
  janitor::clean_names() |> 
  st_transform(crs = st_crs(roads_1))

# Checking the results
# ggplot() +
#   geom_sf(data = region_map) +
#   geom_sf(data = roads_1)
# 
# ggplot() +
#   geom_sf(data = region_map) +
#   geom_sf(data = roads_1 |> st_intersection(region_map))

# Select Roads and Facilties that are within the selected region
roads_1 <- roads_1 |> st_intersection(region_map)
roads_2 <- roads_2 |> st_intersection(region_map)

# Health facilities are as points, polygons and multi-polygons. 
# Lets convert them all into points, since on a large map 
# like a city or state, all will anyways appear as a point

# -------------------------------------------------------------------
# Extracting the Larger Health Facilities in a nicely formatted clean
# sf object, to make plotting easier subsequently--------------------
temp1 <- health_large$osm_points |> 
  select(name, geometry, "addr:district") |> 
  # filter(!is.na(name)) |> 
  mutate(facility_type = "medium")

temp2 <- health_large$osm_polygons |> 
  st_centroid() |> 
  select(name, geometry, "addr:district") |> 
  # filter(!is.na(name)) |> 
  mutate(facility_type = "large")

temp3 <- health_large$osm_multipolygons |> 
  st_centroid() |> 
  select(name, geometry, "addr:district") |> 
  # filter(!is.na(name)) |> 
  mutate(facility_type = "large")

health_large <- bind_rows(
  temp1,
  temp2,
  temp3
)
# Clean the workspace
rm(temp1, temp2, temp3)

# Retain only the health facilities that are within our region
health_large <- health_large |> 
  st_intersection(region_map)

# Now the MAGIC of {sf}: mark the facilities that are near highways
health_large <- health_large |> 
  mutate(
    near_highway = health_large |> 
      
      # Checking which points are within a specific distance of highway
      # The awesome power & speed of sf::st_is_within_distance() 
      st_is_within_distance(
        y = roads_1,
        dist = 200,
        sparse = F
      ) |> 
      as_tibble() |> 
      # Keep value TRUE if facility is near ANY highway
      mutate(
        any_true = rowSums(across(everything(), as.logical)) > 0
      ) |> 
      pull(any_true)
  ) |> 
  mutate(
    near_highway = if_else(
      near_highway,
      "near_highway",
      "far_from_highway"
    )
  )

# Do the same process for snmaller health facilities ----------------
temp1 <- health_small$osm_points |> 
  select(name, geometry, "addr:district") |> 
  # filter(!is.na(name)) |> 
  mutate(facility_type = "small")

temp2 <- health_small$osm_polygons |> 
  st_centroid() |> 
  select(name, geometry, "addr:district") |> 
  # filter(!is.na(name)) |> 
  mutate(facility_type = "medium")

health_small <- bind_rows(
  temp1,
  temp2
)

rm(temp1, temp2)

health_small <- health_small |> 
  st_intersection(region_map)

health_small <- health_small |> 
  mutate(
    near_highway = health_small |> 
      
      # Checking which points are within a specific distance of highway
      st_is_within_distance(
        y = roads_1,
        dist = 200,
        sparse = F
      ) |> 
      as_tibble() |> 
      mutate(
        any_true = rowSums(across(everything(), as.logical)) > 0
      ) |> 
      pull(any_true)
  ) |> 
  mutate(
    near_highway = if_else(
      near_highway,
      "near_highway",
      "far_from_highway"
    )
  )

# Combining health_large and health_small to map shape to facility size
health_combined <- bind_rows(
  health_large,
  health_small
  ) |> 
  mutate(
    facility_type = fct(
      facility_type,
      levels = c("small", "medium", "large")
    )
  ) |> 
  rename(district = addr.district) |> 
  filter(!is.na(district)) |> 
  mutate(
    facility_type = str_to_title(facility_type),
    near_highway = snakecase::to_title_case(near_highway)
  )


rm(health_large, health_small)

# The plot
g <- ggplot() +
  
  # Background map of districts
  geom_sf(
    data = haryana_districts_map,
    linewidth  = 2,
    colour = "grey90",
    alpha = 0.1,
    fill = "transparent"
  ) +
  
  # Background overall map of Haryana
  geom_sf(
    data = region_map,
    linewidth  = 2,
    colour = "grey10",
    fill = "transparent"
  ) +
  
  # Plotting the Highways and other major roads
  geom_sf(
    data = roads_1, 
    linewidth = 0.9, 
    alpha = 1,
    colour = "#702F8AFF" 
  ) +
  geom_sf(
    data = roads_2,
    linewidth = 0.5,
    alpha = 0.8,
    colour = "#9063CDFF"
  ) +
  
  # Plotting the Health Facilties
  geom_sf(
    data = health_combined,
    mapping = aes(
      colour = near_highway,
      size = facility_type,
      shape = facility_type
    ),
    alpha = 0.5
  ) +
  
  scale_size_manual(
    values = c(14, 8, 4)
  ) +
  scale_shape_manual(
    values = c(0, 2, 1)
  ) +
  scale_colour_manual(
    values = c("#FFA400FF", "#862633FF")
  ) +
  
  # Add North-Arrow and Scale
  ggspatial::annotation_north_arrow(
    style = ggspatial::north_arrow_orienteering(
      line_col = "grey20",
      fill = c("grey90", "grey20"),
      text_size = bts,
      text_family = "body_font"
    ),
    location = "tr",
    height = unit(5, "cm"),
    width = unit(5, "cm")
  ) +
  ggspatial::annotation_scale(
    bar_cols = c("grey30", bg_col),
    location = "bl",
    height = unit(0.75, "cm"),
    text_cex = 10,
    text_family = "body_font"
  ) +
  labs(
    title = plot_title,
    subtitle = str_wrap(plot_subtitle, 70),
    caption = plot_caption,
    shape = NULL,
    size = NULL,
    colour = NULL
  ) +
  guides(
    colour = guide_legend(
      override.aes = list(
        size = 10
      )
    )
  ) +
  ggthemes::theme_map(
    base_size = bts,
    base_family = "body_font"
  ) +
  theme(
    # Overall
    plot.margin = margin(10,10,10,10, "mm"),
    plot.title.position = "plot",
    text = element_text(
      colour = text_col,
      lineheight = 0.3,
      margin = margin(0,0,0,0, "mm"),
      hjust = 0
    ),
    
    # Legend
    legend.position = "inside",
    legend.position.inside = c(0, 0.8),
    legend.direction = "vertical",
    legend.text = element_text(
      margin = margin(5,5,5,5, "mm")
    ),
    legend.margin = margin(0,0,0,0, "mm"),
    legend.box = "horizontal",
    legend.box.just = "top",
    legend.key.spacing = unit(10, "mm"),
    
    # Labels
    plot.title = element_text(
      size = 3 * bts,
      hjust = 0.5,
      colour = text_hil,
      margin = margin(30,0,10,0, "mm"),
      family = "title_font"
    ),
    plot.subtitle = element_text(
      size = 0.8 * bts,
      hjust = 0,
      colour = text_hil,
      margin = margin(10,0,-50,0, "mm"),
      family = "body_font"
    ),
    plot.caption = element_textbox(
      family = "caption_font",
      hjust = 0.5,
      margin = margin(10,0,10,0, "mm"),
      colour = text_hil,
      size = 0.5 * bts
    )
  )

###### Computing Data for Table
plot_data_df <- health_combined |> 
  st_drop_geometry() |> 
  as_tibble() |> 
  filter(!(district %in% c("Chandigarh", "Gautam Buddha Nagar",
                         "Haryana", "", " ",
                         "Sahibzada Ajit Singh Nagar"))) |> 
  mutate(
    district = case_when(
      district == "Gurgaon" ~ "Gurugram",
      district == "Nuh" ~ "Mewat",
      .default = district
    )
  ) |> 
  count(district, facility_type, near_highway) |> 
  group_by(district, facility_type) |> 
  mutate(total = sum(n)) |> 
  ungroup()
  # pivot_wider(
  #   id_cols = c(district, facility_type),
  #   names_from = near_highway,
  #   values_from = n,
  #   values_fill = 0
  # ) |> 
  # mutate(
  #   total = `Far from Highway` + `Near Highway`,
  #   perc_on_hwy = `Near Highway` / total,
  #   district = as.character(district),
  #   .keep = "unused"
  # )

strip_labels <- c(
      "Small: Dentists, Pharmacies",
      "Medium: Hospitals & Clinics",
      "Large Institutions"
    ) |> 
  str_wrap(30)
names(strip_labels) <- c("Small", "Medium", "Large")

# Deciding order of district (preference is for descending order)
order_districts <- plot_data_df |> 
  group_by(district) |> 
  summarise(n = sum(n)) |> 
  arrange(n) |> 
  pull(district)

inset_plot <- plot_data_df |> 
  mutate(district = fct(district, levels = order_districts)) |> 
  ggplot(
    mapping = aes(
      y = district,
      x = n
    )
  ) +
  geom_col(
    mapping = aes(
      fill = near_highway
    ),
    alpha = 0.7,
    position = position_stack()
  ) +
  geom_text(
    mapping = aes(
      label = n,
      colour = near_highway
    ),
    hjust = 1.4,
    family = "caption_font",
    size = bts / 20,
    position = position_stack()
  ) +
  geom_text(
    mapping = aes(
      label = total,
      x = total
    ),
    hjust = -0.2,
    family = "body_font",
    size = bts / 12,
    colour = text_col,
    fontface = "bold"
  ) +
  scale_fill_manual(
    values = c("#FFA400FF", "#862633FF")
  ) +
  scale_colour_manual(
    values = darken(c("#FFA400FF", "#862633FF"), 0.5)
  ) +
  scale_y_discrete(
    expand = expansion(0)
  ) +
  scale_x_continuous(
    expand = expansion(c(0, 0.25))
  ) +
  labs(
    y = NULL,
    x = "Number of Health-Care facilties"
  ) +
  facet_wrap(
    ~ facility_type,
    labeller = labeller(
      facility_type = strip_labels
    ),
    scales = "free_x"
  ) +
  coord_cartesian(clip = "off") +
  theme_minimal(
    base_size = bts / 3,
    base_family = "body_font"
  ) +
  theme(
    legend.position = "none",
    plot.margin = margin(0,0,0,0, "mm"),
    plot.background = element_rect(
      fill = "transparent",
      colour = "transparent"
    ),
    panel.background = element_rect(
      fill = "transparent",
      colour = "transparent"
    ),
    panel.grid = element_blank(),
    panel.grid.major.x = element_line(
      colour = alpha(text_col, 0.2),
      linetype = 1,
      linewidth = 0.3
    ),
    axis.line.x = element_line(
      colour = text_hil,
      linewidth = 0.3
    ),
    axis.text.x = element_text(
      margin = margin(5,0,0,0, "mm")
    ),
    axis.text.y = element_text(
      margin = margin(0,5,0,0, "mm")
    ),
    axis.ticks = element_blank(),
    axis.ticks.length = unit(0, "mm")
  )

# Temporary Code to check the output of inset plot
# ggsave(
#   plot = inset_plot,
#   filename = here("projects", "images", 
#                   "highways_health_temp1.png"),
#   height = 8,
#   width = 8,
#   units = "in",
#   bg = bg_col
# )

# Compose the final plot

g2 <- g +
  inset_element(
    p = inset_plot,
    left = 0, right = 0.4,
    bottom = 0.1 , top = 0.3,
    align_to = "full",
    on_top = TRUE,
    clip = FALSE
  )

ggsave(
  plot = g2,
  filename = here("projects", "images", 
                  "highways_health_haryana.png"),
  height = 36,
  width = 24,
  units = "in",
  bg = bg_col
)

Figure 3: The map of Haryana, and its 22 districts, showing the National Highways, and other state highways. Also shown are health facilities (Large Institutions, Medium-sized hospitals, and small sized clinics and pharmacies. They are colour coded by their proximity to National Highways (defined as a distance of less than 200 metres). The inset bar-graph shows the number of health-facilities of each size in each district.

Part 7: Final Product for India (in-progress)

Code

# Basics of Final Visualization: Fonts
font_add_google("Rouge Script", "title_font")
font_add_google("El Messiri", "body_font")
font_add_google("Saira Extra Condensed", "caption_font")
showtext_auto()

text_col <- "grey10"
text_hil <- "grey20"
bts <- 150
bg_col <- "white"

# Caption stuff for the plot
sysfonts::font_add(
  family = "Font Awesome 6 Brands",
  regular = here::here("docs", "Font Awesome 6 Brands-Regular-400.otf")
)
github <- "&#xf09b"
github_username <- "aditya-dahiya"
xtwitter <- "&#xe61b"
xtwitter_username <- "@adityadahiyaias"
social_caption_1 <- glue::glue("<span style='font-family:\"Font Awesome 6 Brands\";'>{github};</span> <span style='color:{text_hil}'>{github_username}  </span>")
social_caption_2 <- glue::glue("<span style='font-family:\"Font Awesome 6 Brands\";'>{xtwitter};</span> <span style='color:{text_hil}'>{xtwitter_username}</span>")

plot_title <- "India: Health-Care Institutions"

plot_caption <- paste0(
  "**Data:** Open Street Maps; Survey of India", 
  " |  **Code:** ", 
  social_caption_1, 
  " |  **Graphics:** ", 
  social_caption_2
  )

rm(github, github_username, xtwitter, 
   xtwitter_username, social_caption_1, 
   social_caption_2)

# Get pre-saved data (Trial Phase - use pre-saved data to 
# save laoding times)
# Note: I have split downloading health facilties into larger and
#       Smaller ones to save on downloading times, and prevent 
#       breakage during downloading such large datasets
health_small <- readRDS(file = "india_health_small.rds")
health_large <- readRDS(file = "india_health_large.rds")
# roads_1 <- readRDS(file = "india_roads_1.rds")
# roads_2 <- readRDS(file = "haryana_roads_2.rds")
# 
# # Get only the relevant sf object out of each downlaoded data-set
# # For roads, only "LINES" are relevant sf feature.
# roads_1 <- roads_1$osm_lines
# roads_2 <- roads_2$osm_lines

# Get a Map of the State / City and keep only the highways and 
# facilties within that State / City for plotting (instead of entire)
# bounding box
region_map <- sf::read_sf(
  here::here("data", "india_map", "India_Country_Boundary.shp")
) |> 
  janitor::clean_names() |> 
  st_transform(crs = st_crs(health_large$osm_points))

# # Select Roads and Facilties that are within the selected region
# roads_1 <- roads_1 |> st_intersection(region_map)
# roads_2 <- roads_2 |> st_intersection(region_map)

# Health facilities are as points, polygons and multi-polygons. 
# Lets convert them all into points, since on a large map 
# like a city or state, all will anyways appear as a point

# ----------------------------------------------------------
# Extracting the Larger Health Facilities in a nicely formatted clean
# sf object, to make plotting easier subsequently--------------------
temp1 <- health_large$osm_points |> 
  select(name, geometry, "addr:state") |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "medium")

temp2 <- health_large$osm_polygons |> 
  st_centroid() |> 
  select(name, geometry, "addr:state") |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "large")

sf_use_s2(FALSE)
temp3 <- health_large$osm_multipolygons |> 
  st_centroid(of_largest_polygon = TRUE) |> 
  select(name, geometry, "addr:state") |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "large")

health_large <- bind_rows(
  temp1,
  temp2,
  temp3
)
# Clean the workspace
rm(temp1, temp2, temp3)
sf_use_s2(TRUE)
# Retain only the health facilities that are within our region
health_large <- health_large |> 
  st_intersection(region_map)

# Now the MAGIC of {sf}: mark the facilities that are near highways
# health_large <- health_large |> 
#   mutate(
#     near_highway = health_large |> 
#       
#       # Checking which points are within a specific distance of highway
#       # The awesome power & speed of sf::st_is_within_distance() 
#       st_is_within_distance(
#         y = roads_1,
#         dist = 200,
#         sparse = F
#       ) |> 
#       as_tibble() |> 
#       # Keep value TRUE if facility is near ANY highway
#       mutate(
#         any_true = rowSums(across(everything(), as.logical)) > 0
#       ) |> 
#       pull(any_true)
#   ) |> 
#   mutate(
#     near_highway = if_else(
#       near_highway,
#       "near_highway",
#       "far_from_highway"
#     )
#   )

# Do the same process for snmaller health facilities ----------------
temp1 <- health_small$osm_points |> 
  select(name, geometry, "addr:state") |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "small")

temp2 <- health_small$osm_polygons |> 
  st_centroid() |> 
  select(name, geometry, "addr:state") |> 
  filter(!is.na(name)) |> 
  mutate(facility_type = "medium")

health_small <- bind_rows(
  temp1,
  temp2
)

rm(temp1, temp2)

health_small <- health_small |> 
  st_intersection(region_map)

# health_small <- health_small |> 
#   mutate(
#     near_highway = health_small |> 
#       
#       # Checking which points are within a specific distance of highway
#       st_is_within_distance(
#         y = roads_1,
#         dist = 200,
#         sparse = F
#       ) |> 
#       as_tibble() |> 
#       mutate(
#         any_true = rowSums(across(everything(), as.logical)) > 0
#       ) |> 
#       pull(any_true)
#   ) |> 
#   mutate(
#     near_highway = if_else(
#       near_highway,
#       "near_highway",
#       "far_from_highway"
#     )
#   )

# Combining health_large and health_small to map shape to facility size
health_combined <- bind_rows(
  health_large,
  health_small
  ) |> 
  mutate(
    facility_type = fct(
      facility_type,
      levels = c("small", "medium", "large")
    )
  ) |> 
  rename(district = addr.district) |> 
  filter(!is.na(district)) |> 
  mutate(
    facility_type = str_to_title(facility_type),
    near_highway = snakecase::to_title_case(near_highway)
  )


rm(health_large, health_small)

# The plot
g <- ggplot() +
  
  # Background map of districts
  geom_sf(
    data = haryana_districts_map,
    linewidth  = 2,
    colour = "grey90",
    alpha = 0.1,
    fill = "transparent"
  ) +
  
  # Background overall map of Haryana
  geom_sf(
    data = region_map,
    linewidth  = 2,
    colour = "grey10",
    fill = "transparent"
  ) +
  
  # Plotting the Highways and other major roads
  geom_sf(
    data = roads_1, 
    linewidth = 0.9, 
    alpha = 0.8,
    colour = "#702F8AFF" 
  ) +
  geom_sf(
    data = roads_2,
    linewidth = 0.5,
    alpha = 0.5,
    colour = "#9063CDFF"
  ) +
  
  # Plotting the Health Facilties
  geom_sf(
    data = health_combined,
    mapping = aes(
      colour = near_highway,
      size = facility_type,
      shape = facility_type
    ),
    alpha = 0.5
  ) +
  
  scale_size_manual(
    values = c(10, 7, 3)
  ) +
  scale_shape_manual(
    values = c(0, 2, 20)
  ) +
  scale_colour_manual(
    values = c("#FFA400FF", "#862633FF")
  ) +
  
  # Add North-Arrow and Scale
  ggspatial::annotation_north_arrow(
    style = ggspatial::north_arrow_orienteering(
      line_col = "grey20",
      fill = c("grey90", "grey20"),
      text_size = bts,
      text_family = "body_font"
    ),
    location = "tr",
    height = unit(5, "cm"),
    width = unit(5, "cm")
  ) +
  ggspatial::annotation_scale(
    bar_cols = c("grey30", bg_col),
    location = "bl",
    height = unit(0.75, "cm"),
    text_cex = 10,
    text_family = "body_font"
  ) +
  labs(
    title = plot_title,
    caption = plot_caption,
    shape = NULL,
    size = NULL,
    colour = NULL
  ) +
  guides(
    colour = guide_legend(
      override.aes = list(
        size = 10
      )
    )
  ) +
  ggthemes::theme_map(
    base_size = bts,
    base_family = "body_font"
  ) +
  theme(
    # Overall
    plot.margin = margin(10,10,10,10, "mm"),
    plot.title.position = "plot",
    text = element_text(
      colour = text_col,
      lineheight = 0.3,
      margin = margin(0,0,0,0, "mm"),
      hjust = 0
    ),
    
    # Legend
    legend.position = "inside",
    legend.position.inside = c(0, 0.8),
    legend.direction = "vertical",
    legend.text = element_text(
      margin = margin(5,5,5,5, "mm")
    ),
    legend.margin = margin(0,0,0,0, "mm"),
    legend.box = "horizontal",
    legend.box.just = "top",
    legend.key.spacing = unit(10, "mm"),
    
    # Labels
    plot.title = element_text(
      size = 3 * bts,
      hjust = 0.5,
      colour = text_hil,
      margin = margin(30,0,30,0, "mm"),
      family = "title_font"
    ),
    plot.caption = element_textbox(
      family = "caption_font",
      hjust = 0.5,
      margin = margin(10,0,10,0, "mm"),
      colour = text_hil,
      size = 0.75 * bts
    )
  )

###### Computing Data for Table
plot_data_df <- health_combined |> 
  st_drop_geometry() |> 
  as_tibble() |> 
  count(district, facility_type, near_highway) |> 
  pivot_wider(
    id_cols = c(district, facility_type),
    names_from = near_highway,
    values_from = n,
    values_fill = 0
  ) |> 
  mutate(
    total = `Far from Highway` + `Near Highway`,
    perc_on_hwy = `Near Highway` / total,
    district = as.character(district),
    .keep = "unused"
  ) |> 
  filter(!(district %in% c("Chandigarh", "Gautam Buddha Nagar",
                         "Haryana", "", " ",
                         "Sahibzada Ajit Singh Nagar"))) |> 
  mutate(
    district = case_when(
      district == "Gurgaon" ~ "Gurugram",
      district == "Nuh" ~ "Mewat",
      .default = district
    )
  )

strip_labels <- c(
      "Small facilties: Dentists, Pharmacies",
      "Medium: Clinics, Doctors's Practices",
      "Large: Hospitals, Nursing Homes, Institutions"
    ) |> 
  str_wrap(30)
names(strip_labels) <- c("Small", "Medium", "Large")

inset_plot <- plot_data_df |> 
  ggplot(
    mapping = aes(
      y = district,
      x = total
    )
  ) +
  geom_col(
    alpha = 0.5
  ) +
  geom_text(
    mapping = aes(
      label = paste0("(", round(perc_on_hwy * 100, 2), 
                     "% )"),
      x = 0
    ),
    nudge_x = 0.01,
    hjust = 0
  ) +
  geom_text(
    mapping = aes(
      label = paste0(total)
    ),
    nudge_x = 0.005,
    hjust = 0
  ) +
  labs(
    y = NULL,
    x = "Health facilities located on the Highway"
  ) +
  facet_wrap(
    ~ facility_type,
    labeller = labeller(
      facility_type = strip_labels
    ),
    scales = "free_x"
  ) +
  scale_x_continuous(
    expand = expansion(c(0, 0.25))
  ) +
  theme_minimal(
    base_size = 10
  ) +
  theme()


inset_plot

ggsave(
  plot = g,
  filename = here("projects", "images", 
                  "highways_health_haryana.png"),
  height = 36,
  width = 24,
  units = "in",
  bg = bg_col
)

Part 8 : Thumbnail for webpage, and list of packages used

Code

# Saving a thumbnail

library(magick)
# Saving a thumbnail for the webpage
image_read(here::here("projects", "images", 
                      "highways_health_haryana.png")) |> 
  image_resize(geometry = "400") |> 
  image_write(
    here::here(
      "projects", 
      "images", 
      "highways_health_haryana_thumbnail.png"
    )
  )

Code

# Data Import and Wrangling Tools
library(tidyverse)            # All things tidy
library(here)                 # File locations and paths

# Final plot tools
library(scales)               # Nice Scales for ggplot2
library(fontawesome)          # Icons display in ggplot2
library(ggtext)               # Markdown text support for ggplot2
library(showtext)             # Display fonts in ggplot2
library(colorspace)           # Lighten and Darken colours
library(patchwork)            # Combining plots
library(gt)                   # Beautiful Tables in R

# Mapping tools
library(sf)                   # All spatial objects in R
library(osmdata)              # Getting Open Street Maps data

sessioninfo::session_info()$packages |> 
  as_tibble() |> 
  select(package, 
         version = loadedversion, 
         date, source) |> 
  arrange(package) |> 
  janitor::clean_names(
    case = "title"
  ) |> 
  gt::gt() |> 
  gt::opt_interactive(
    use_search = TRUE
  )

Table 1: R Packages and their versions used in the creation of this page and graphics

References

Mark Padgham, Bob Rudis, Robin Lovelace, and Maëlle Salmon. 2017. “Osmdata” 2: 305. https://doi.org/10.21105/joss.00305.

Pebesma, Edzer, and Roger Bivand. 2023. “Spatial Data Science: With Applications in r.” https://doi.org/10.1201/9780429459016.