Customizing {geofacet} for plotting geographically faceted graphs

Exploring how faceted plots cn be re-arragned to roughly follow the geographical location of sub-units.

Part 1: Working with the {geofacet} package

The geofacet package provides an innovative approach to visualizing geographically structured data by arranging small multiples (facets) in a grid that mirrors the spatial arrangement of regions on a map. This facilitates comparison across regions while retaining geographical context. The package includes built-in grids for various regions, such as state_grid1 for the United States, but users can also create custom grids using grid_preview() to visualize the structure. Integration with ggplot2 allows for seamless creation of geofaceted visualizations, making it possible to apply any ggplot-compatible data. The package is highly customizable, with functions like facet_geo() providing flexibility to specify grids and scales. Additionally, geofacet supports data from other countries, such as Japan and Australia, with pre-built grids available for quick exploration. For more advanced applications, users can define their own layouts using add_custom_grid(). This makes geofacet a powerful tool for researchers and analysts working with spatially contextual data.

Install package, load required libraries

# install.packages("geofacet")
library(tidyverse)             # Data Wrangling
library(sf)                    # Simple Features in R
library(geofacet)              # Geographic faceting
library(patchwork)             # Compiling Plots
library(gt)                    # Great Tables in R  


# Setting themes etc. ----------------------------------
sysfonts::font_add_google("Saira Condensed", "body_font")
showtext::showtext_auto()
theme_set(
  theme_minimal(
    base_family = "body_font",
    base_size = 14,
    base_line_size = 0.2
  ) +
    theme(
      text = element_text(
        colour = "grey20"
      ),
      plot.title = element_text(
        size = 24,
        margin = margin(0.1, 0.1, 0.1, 0.1, "lines")
      ),
      plot.subtitle = element_text(
        size = 16,
        margin = margin(0.1, 0.1, 0.1, 0.1, "lines"),
        lineheight = 0.3
      )
    )
)

Exploring the package with an example

This code demonstrates how to use the geofacet package to create geofaceted bar plots shown in Figure 1, of state ranks for U.S. states, utilizing two different in-built grid layouts: us_state_grid1 and us_state_grid2. It showcases comparison across layouts by combining plots with patchwork, adding annotations, and saving the final visualization as an image.

data("state_ranks")

g1 <- state_ranks |> 
  as_tibble() |> 
  mutate(variable = str_to_title(variable)) |> 
  ggplot(aes(rank, variable, fill = variable)) +
  geom_col() +
  geom_text(
    aes(label = rank),
    size = 2,
    colour = "white",
    hjust = 1.1
  ) +
  facet_geo(~state) +
  paletteer::scale_fill_paletteer_d(
    "ggthemes::Classic_Purple_Gray_6"
  ) +
  labs(
    y = NULL, x = "Rank of the State",
    title = "Layout 1: `us_state_grid1`"
  ) +
  theme_minimal(
    base_family = "body_font"
  ) +
  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    plot.title = element_text(
      hjust = 0, size = 48
    )
  )

g2 <- state_ranks |> 
  as_tibble() |> 
  mutate(variable = str_to_title(variable)) |> 
  ggplot(aes(rank, variable, fill = variable)) +
  geom_col() +
  geom_text(
    aes(label = rank),
    size = 2,
    colour = "white",
    hjust = 1.1
  ) +
  facet_geo(~state, grid = "us_state_grid2") +
  paletteer::scale_fill_paletteer_d(
    "ggthemes::Classic_Purple_Gray_6"
  ) +
  labs(
    y = NULL, x = "Rank of the State",
    title = "Layout 2: `us_state_grid2`"
  ) +
  theme_minimal(
    base_family = "body_font"
  ) +
  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    axis.text.x = element_blank(),
    plot.title = element_text(
      hjust = 0, size = 48
    )
  )

# Compiling Plots ----------------------------------------
g <- g1 + g2 + 
  plot_layout(
    nrow = 1,
    guides = "collect"
  )  +
  plot_annotation(
    title = "Two in-built payouts from {geofacet}",
    subtitle = "Displaying state ranks for US states, using in-built dataset from {geofacet}",
    tag_levels = "a",
    tag_prefix = "(",
    tag_suffix = ")",
    theme = theme(
      plot.title = element_text(
        hjust = 0.5, size = 54,
        margin = margin(2,0,5,0, "pt"),
        face = "bold"
      ),
      plot.subtitle = element_text(
        hjust = 0.5, size = 36,
        margin = margin(2,0,5,0, "pt")
      )
    )
  ) &
  theme(
    plot.tag = element_text(
      size = 48, face = "bold"
    ),
    axis.ticks = element_blank(),
    axis.ticks.length = unit(0, "pt"),
    axis.text = element_text(
      margin = margin(0,0,0,0, "mm")
    ),
    panel.spacing = unit(0, "pt"),
    strip.text = element_text(
      margin = margin(0,0,0,0, "pt"),
      size = 24, face = "bold"
    )
  )

ggsave(
  plot = g,
  filename = here::here("geocomputation", "images", "custom_geofacet_1.png"),
  height = 1300,
  width = 3000,
  unit = "px",
  bg = "white"
)

Figure 1: Maps of USA with statitics (ranks) for each state plotted as horizontal bar charts. Facet layouts compiled using: for map (a) geofacet::facet_geo() and for map (b) using geofacet::facet_geo(grid = "us_state_grid2")

Part 2: An Example submitting a custom Grid to {geofacet}

About the 22 districts in the State of Haryana

This code creates a detailed summary table of the 22 districts of Haryana, including their names, codes, and approximate centroid locations (latitude and longitude). Using the tibble package, it structures the data, and the gt package formats it into an interactive table with customized column labels and number formatting, shown in Table 1.

haryana_districts <- tibble(
  District_Name = c(
    "Ambala", "Bhiwani", "Charkhi Dadri", "Faridabad", "Fatehabad",
    "Gurugram", "Hisar", "Jhajjar", "Jind", "Kaithal",
    "Karnal", "Kurukshetra", "Mahendragarh", "Nuh", "Palwal",
    "Panchkula", "Panipat", "Rewari", "Rohtak", "Sirsa",
    "Sonipat", "Yamunanagar"
  ),
  Two_Letter_Code = c(
    "AM", "BH", "CD", "FR", "FT", "GU", "HI", "JH", "JI", "KT",
    "KR", "KU", "MA", "NU", "PW", "PK", "PP", "RE", "RO", "SI",
    "SN", "YN"
  ),
  Three_Letter_Code = c(
    "AMB", "BHW", "CHD", "FAR", "FAT", "GUR", "HIS", 
    "JHA", "JIN", "KAI",
    "KAR", "KUR", "MAH", "NUH", "PAL", "PAN", "PNP", 
    "REW", "ROH", "SIR", "SON", "YAM"
  ),
  Latitude = c(
    30.3782, 28.7930, 28.5921, 28.4089, 29.5252,
    28.4595, 29.1492, 28.6064, 29.3162, 29.8010,
    29.6857, 29.9695, 28.2692, 28.1070, 28.1447,
    30.6942, 29.3909, 28.1970, 28.8955, 29.5349,
    28.9931, 30.1290
  ),
  Longitude = c(
    76.7767, 76.1390, 76.2711, 77.3178, 75.4540,
    77.0266, 75.7217, 76.6565, 76.3144, 76.3995,
    76.9905, 76.8783, 76.1521, 77.0010, 77.3260,
    76.8606, 76.9635, 76.6170, 76.6066, 75.0280,
    77.0151, 77.2674
  )
)

gt(haryana_districts) |> 
  cols_label_with(fn = snakecase::to_title_case) |> 
  fmt_number(columns = c(Latitude, Longitude)) |> 
  opt_interactive() |> 
  gtExtras::gt_theme_538()

Table 1: Summary of 22 Districts in the state of Haryana, along with their approximate centroid locations

Creating Haryana’s Custom {geofacet} grid

This code creates a custom geofacet grid for visualizing Haryana’s districts in a spatially relevant layout. It defines the grid’s structure using a data frame, previews the layout with grid_preview(), and customizes the appearance with titles and themes. Finally, the grid is saved as an image, shown in Figure 2.

haryana_grid <- data.frame(
  name = c("Panchkula", "Ambala", "Kaithal", "Yamunanagar", "Fatehabad", "Kurukshetra", "Karnal", "Sirsa", "Jind", "Panipat", "Hisar", "Rohtak", "Bhiwani", "Sonipat", "Charkhi Dadri", "Jhajjar", "Rewari", "Mahendragarh", "Faridabad", "Gurugram", "Nuh", "Palwal"),
  code = c("PAN", "AMB", "KAI", "YAM", "FAT", "KUR", "KAR", "SIR", "JIN", "PNP", "HIS", "ROH", "BHW", "SON", "CHD", "JHA", "REW", "MAH", "FAR", "GUR", "NUH", "PAL"),
  row = c(1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 6, 6, 7, 7, 7, 8, 8),
  col = c(4, 4, 3, 5, 2, 3, 4, 1, 3, 4, 2, 4, 2, 5, 3, 4, 3, 2, 5, 4, 4, 5),
  stringsAsFactors = FALSE
)

g <- geofacet::grid_preview(
  haryana_grid,
  label = "name"
  ) +
  labs(
    title = "A custom geofacet grid for the districts in State of Haryana (India)",
    subtitle = "Created using geofacet::grid_design()",
    x = "Column Number", y = "Row Number"
  ) +
  theme(
    plot.title = element_text(
      hjust = 0.5
    ),
    plot.subtitle = element_text(
      hjust = 0.5
    )
  )

ggsave(
  plot = g,
  filename = here::here("geocomputation", "images", 
                        "custom_geofacet_2.png"),
  height = 700,
  width = 800,
  unit = "px",
  bg = "white"
)

Figure 2: A custom grid created using geofacet::design_grid() is displayed using geofacet::grid_preview()

An example dataset to plot

This code creates an interactive table summarizing the percentage area of each zone, categorized by proximity to health-care facilities, across districts in Haryana. For how this data was calculated, and the compete project, check out this webpage. It uses tidyverse for data manipulation, reshapes the dataset with pivot_wider(), and formats the table using gt() with an interactive display and customized headers, shown in Table 2 .

library(tidyverse)
library(gt)
hy_zones <- structure(list(buf_dist = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 
7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 
3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 
6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 
2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 
4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 
1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 
5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 
2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 
4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L), levels = c("< 1 km", 
"1 - 4 km", "4 - 8 km", "8 - 12 km", "12 - 16 km", "16 - 20 km", 
"> 20 km"), class = "factor"), district = c("Ambala", "Ambala", 
"Ambala", "Ambala", "Ambala", "Ambala", "Ambala", "Bhiwani", 
"Bhiwani", "Bhiwani", "Bhiwani", "Bhiwani", "Bhiwani", "Bhiwani", 
"Faridabad", "Faridabad", "Faridabad", "Faridabad", "Faridabad", 
"Faridabad", "Fatehabad", "Fatehabad", "Fatehabad", "Fatehabad", 
"Fatehabad", "Fatehabad", "Fatehabad", "Gurugram", "Gurugram", 
"Gurugram", "Gurugram", "Gurugram", "Gurugram", "Hisar", "Hisar", 
"Hisar", "Hisar", "Hisar", "Hisar", "Hisar", "Jhajjar", "Jhajjar", 
"Jhajjar", "Jhajjar", "Jhajjar", "Jhajjar", "Jind", "Jind", "Jind", 
"Jind", "Jind", "Jind", "Jind", "Kaithal", "Kaithal", "Kaithal", 
"Kaithal", "Kaithal", "Kaithal", "Kaithal", "Karnal", "Karnal", 
"Karnal", "Karnal", "Karnal", "Karnal", "Karnal", "Kurukshetra", 
"Kurukshetra", "Kurukshetra", "Kurukshetra", "Kurukshetra", "Kurukshetra", 
"Mahendragarh", "Mahendragarh", "Mahendragarh", "Mahendragarh", 
"Mahendragarh", "Mahendragarh", "Mahendragarh", "Mewat", "Mewat", 
"Mewat", "Mewat", "Mewat", "Mewat", "Palwal", "Palwal", "Palwal", 
"Palwal", "Palwal", "Palwal", "Panchkula", "Panchkula", "Panchkula", 
"Panchkula", "Panchkula", "Panchkula", "Panipat", "Panipat", 
"Panipat", "Panipat", "Panipat", "Panipat", "Rewari", "Rewari", 
"Rewari", "Rewari", "Rewari", "Rewari", "Rohtak", "Rohtak", "Rohtak", 
"Rohtak", "Rohtak", "Rohtak", "Rohtak", "Sirsa", "Sirsa", "Sirsa", 
"Sirsa", "Sirsa", "Sirsa", "Sirsa", "Sonipat", "Sonipat", "Sonipat", 
"Sonipat", "Sonipat", "Sonipat", "Sonipat", "Yamunanagar", "Yamunanagar", 
"Yamunanagar", "Yamunanagar", "Yamunanagar", "Yamunanagar", "Yamunanagar", 
"Charkhi Dadri", "Charkhi Dadri", "Charkhi Dadri", "Charkhi Dadri", 
"Charkhi Dadri", "Charkhi Dadri"), perc = c(0.0329, 0.1599, 0.2706, 
0.2768, 0.2065, 0.0467, 0.0066, 0.0138, 0.1099, 0.2747, 0.3029, 
0.2181, 0.0781, 0.0025, 0.1514, 0.3398, 0.3459, 0.1553, 0.0076, 
0, 0.0131, 0.1151, 0.2478, 0.3059, 0.2378, 0.0765, 0.0039, 0.134, 
0.3348, 0.3725, 0.144, 0.0148, 0, 0.0198, 0.1407, 0.3248, 0.3301, 
0.1507, 0.0327, 0.0011, 0.0263, 0.2238, 0.4302, 0.2664, 0.0528, 
5e-04, 0.0152, 0.1458, 0.2934, 0.2476, 0.1642, 0.0774, 0.0564, 
0.0164, 0.1475, 0.3411, 0.2832, 0.1225, 0.0676, 0.0218, 0.0247, 
0.2134, 0.4162, 0.219, 0.0831, 0.0358, 0.0077, 0.0358, 0.2324, 
0.4168, 0.2253, 0.0823, 0.0073, 0.0226, 0.167, 0.3579, 0.3207, 
0.0972, 0.0212, 0.0134, 0.0167, 0.1676, 0.4231, 0.2727, 0.1134, 
0.0065, 0.0242, 0.1469, 0.3006, 0.3515, 0.1731, 0.0037, 0.0618, 
0.266, 0.3477, 0.2475, 0.0607, 0.0162, 0.0445, 0.2475, 0.4688, 
0.2171, 0.0217, 4e-04, 0.0262, 0.2104, 0.4112, 0.269, 0.0819, 
0.0013, 0.0359, 0.1939, 0.3301, 0.2855, 0.1055, 0.0407, 0.0082, 
0.0088, 0.088, 0.2127, 0.2543, 0.2169, 0.1081, 0.1112, 0.0313, 
0.2556, 0.4207, 0.2006, 0.0532, 0.0293, 0.0092, 0.029, 0.1784, 
0.3443, 0.2619, 0.1442, 0.0403, 0.0019, 0.0246, 0.2249, 0.4622, 
0.2321, 0.0512, 0.005)), row.names = c(NA, -144L), class = c("tbl_df", 
"tbl", "data.frame")) |> 
  mutate(district = if_else(district == "Mewat", "Nuh", district))

hy_zones |> 
  mutate(perc = perc * 100) |> 
  pivot_wider(
    id_cols = district,
    names_from = buf_dist,
    values_from = perc
  ) |> 
  gt() |> 
  cols_label(
    district = "District"
  ) |> 
  opt_interactive() |> 
  tab_header(
    title = "Percentage area of each district, in each zone",
    subtitle = "Column headings represent the distance from the nearest health-care facility"
  ) |> 
  gtExtras::gt_theme_538()

Table 2: The Raw data for distance zones from nearest health-care facility for the 22 districts of Haryana

Percentage area of each district, in each zone

Column headings represent the distance from the nearest health-care facility

A customized {geofacet} graphic

This code demonstrates how to visualize the percentage of Haryana’s district areas within certain distance zones from health-care facilities using both standard facet_wrap() and customized facet_geo() layouts. The plots use stacked polar bar charts for comparison, with a custom geofacet layout reflecting Haryana’s spatial arrangement. The combined plots are styled, annotated, and saved as a PNG file for further use, shown in Figure 3 .

# Improve hy_zones tibble into final plotting tibble "plotdf"
plotdf <- hy_zones |> 
  left_join(
    haryana_districts |> 
      select(District_Name, Three_Letter_Code),
    by = join_by(district == District_Name)
  ) |> 
  rename(code = Three_Letter_Code,
         name = district)

# Routine facet plot ----------------------------------
g1 <- plotdf |> 
  ggplot(
    mapping = aes(
      x = 1,
      y = perc,
      group = buf_dist,
      fill = buf_dist
    )
  ) +
  geom_col(
    position = position_stack(),
    colour = "white",
    linewidth = 0.2
  ) +
  geom_text(
    mapping = aes(
      label = paste0(round(100 * perc, 1), "%")
    ),
    position = position_stack(
      vjust = 0.5
    ),
    family = "body_font",
    size = 4,
    check_overlap = TRUE
  ) +
  scale_x_continuous(expand = expansion(0)) +
  scale_fill_manual(
    values = paletteer::paletteer_d("MexBrewer::Taurus1")
  ) +
  guides(
    fill = guide_legend(
      nrow = 1
    )
  ) +
  facet_wrap(~name) +
  coord_polar(
    theta = "y"
  ) +
  labs(
    x = NULL, y = NULL,
    fill = "Distance from nearest health-care facility",
    title = "Regular facet_wrap() layout"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_size = 24
  )

# A custom geofacet layout ------------------------------
g2 <- plotdf |> 
  drop_na() |> 
  ggplot(
    mapping = aes(
      x = 1,
      y = perc,
      group = buf_dist,
      fill = buf_dist
    )
  ) +
  geom_col(
    position = position_stack(),
    colour = "white",
    linewidth = 0.2
  ) +
  geom_text(
    mapping = aes(
      label = paste0(round(100 * perc, 1), "%")
    ),
    position = position_stack(
      vjust = 0.5
    ),
    family = "body_font",
    size = 4,
    check_overlap = TRUE
  ) +
  scale_x_continuous(expand = expansion(0)) +
  scale_fill_manual(
    values = paletteer::paletteer_d("MexBrewer::Taurus1")
  ) +
  guides(
    fill = guide_legend(
      nrow = 1
    )
  ) +
  facet_geo(
    ~code, 
    grid = haryana_grid,
    label = "name"
  ) +
  coord_radial(
    theta = "y",
    expand = FALSE
  ) +
  labs(
    x = NULL, y = NULL,
    fill = "Distance from nearest health-care facility",
    title = "Customized facet_geo() layout"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_size = 24
  )


# Compiling Plots ----------------------------------------
g <- g1 + g2 + 
  plot_layout(
    nrow = 1,
    guides = "collect"
  )  +
  plot_annotation(
    title = "Comparing facet_wrap() and facet_geo()",
    subtitle = "Percentage area of each district in Haryana that falls within a certain distance zone from nearest health-care facility",
    tag_levels = "a",
    tag_prefix = "(",
    tag_suffix = ")",
    theme = theme(
      plot.subtitle = element_text(
        hjust = 0.5, size = 36,
        margin = margin(2,0,5,0, "pt")
      )
    )
  ) &
  theme(
    panel.grid = element_blank(),
    axis.text = element_blank(),
    legend.position = "bottom",
    legend.title.position = "top",
    panel.spacing = unit(0, "pt"),
    strip.text = element_text(
      margin = margin(0,0,0,0, "pt")
    ),
    legend.title = element_text(hjust = 0.5, size = 40,
                                margin = margin(0,0,2,0, "pt")),
    legend.text = element_text(
      margin = margin(0,4,0,2, "pt")
    ),
    plot.title.position = "plot",
    plot.tag = element_text(
      size = 48,
      margin = margin(2,2,2,2, "pt")
    ),
    plot.margin = margin(0,5,0,5, "pt"),
    legend.margin = margin(0,0,0,0, "pt"),
    legend.box.margin = margin(0,0,0,0, "pt"),
    plot.title = element_text(
      size = 60, 
      margin = margin(2,0,2,0, "pt"),
      hjust = 0.5
    )
  )

ggsave(
  plot = g,
  filename = here::here("geocomputation", "images", "custom_geofacet_4.png"),
  height = 2400,
  width = 3000,
  unit = "px",
  bg = "white"
)

Figure 3: Two plots. (a) Displaying the percentage area of each district that falls within a certain zone as pie charts arranged in regular faceted pattern of {ggplot2} using facet_wrap(). (b) Same pie charts now arranged in a custom {geofacet} layout provided using facet_geo()

The Figure 4 shows a graphic for custom facet_geo() to publish online. The full code to generate it is here.

Figure 4: Pie-charts depicting the percentage area of each district of Haryana that falls within a particular distance zone from the nearest health-care facility. The charts are arranged in a customized facet pattern using geofacet::facet_geo() to mimic the approximate geographic location of the districts.