Combining with Population Density Rasters with OSM data in R with {sf}, {osmextract} and {ggplot2}

Combining with Population Density Raster Dataset with OSM maps

The population raster data is sourced from Global High-Resolution Annual Population Grids (2000-2023), v1 by Ciesin, Center for International Earth Science Information Network, Columbia University (CIESIN, 2024). This dataset provides high-resolution population estimates globally, integrating census and administrative data with geospatial modeling to refine population distributions. The data, maintained by NASA’s Socioeconomic Data and Applications Center (SEDAC), is valuable for demographic analysis, urban planning, and environmental studies. The full dataset and documentation are available on Zenodo.

The Haryana population density map is generated using high-resolution raster data from the Global High-Resolution Annual Population Grids (2000-2023) dataset. The {terra} package is used to handle raster operations, where the rast() function reads the population density raster, followed by crop() and mask() to limit the data to the Haryana state boundary. To ensure effective visualization, values below or equal to zero are replaced with 0.01 for smooth log transformation. The road network is overlaid using {sf}, with geom_sf() displaying highways extracted using {osmextract}. The population raster is plotted using geom_spatraster() from {ggplot2} and styled with a log-transformed color scale via {paletteer} to highlight variations in density.

Comparing with some districts

The inset district maps focus on seven key districts, emphasizing transit-oriented development along highways. The plot_district() function extracts population density data for a given district by using crop() and mask() from {terra}. To enhance clarity, districts with multiple polygons (such as Faridabad) are filtered using st_cast("POLYGON"). The district’s road network is extracted with st_intersection(), ensuring only relevant highways are displayed. {patchwork} is used to arrange the Haryana map alongside its inset districts with a custom layout via plot_layout(), effectively demonstrating the correlation between population density and road infrastructure.

# 2022 year Global Population Density 30 sec arc resolution
# url <- "https://zenodo.org/records/11179644/files/GlobPOP_Count_30arc_2022_I32.tiff?download=1"
# 
output_file <- "GlobPOP_Count_30arc_2022_I32.tiff"
# download.file(url, output_file, mode = "wb")

haryana_pop_rast <- rast(output_file) |> 
  terra::crop(haryana_boundary) |> 
  terra::mask(haryana_boundary)

# Ensure all negative and zero values are replaced with 0.01
# (For easy plotting with log transformation scale)
haryana_pop_rast[haryana_pop_rast <= 0] <- 0.01

g <- ggplot() +
  
  # Population Density Raster
  geom_spatraster(data = haryana_pop_rast) +
  paletteer::scale_fill_paletteer_c(
    "grDevices::YlOrBr",
    direction = -1,
    na.value = "transparent",
    transform = "log",
    limits = c(10, 1e5),
    oob = scales::squish,
    breaks = c(0, 10, 100, 1000, 1e4),
    labels = scales::label_number(big.mark = ",")
  ) +

  # Road Network
  geom_sf(
    data = df1 |> filter(highway %in% c(wid1, wid2)),
    linewidth = 0.3,
    alpha = 0.9
  ) +

  geom_sf(
    data = haryana_boundary,
    linewidth = 1.2,
    colour = "black", 
    fill = NA
  ) +
  labs(
    # title = "Population Density vs. Road Network (Haryana)",
    # subtitle = "Bulk data download using {osmextract}, geocomputation\nwith {sf} and plotting with {ggplot2}",
    # caption = plot_caption,
    fill = "Population Density\n(Persons per sq. km.)"
  ) +
  theme_minimal(
    base_family = "body_font",
    base_size = 40
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.02,0.05),
    legend.justification = c(0,0),
    legend.direction = "vertical",
    text = element_text(
      colour = "grey30",
      lineheight = 0.3,
      hjust = 0.5
    ),
    plot.title = element_text(
      margin = margin(10,0,0,0, "pt"),
      hjust = 0.5,
      size = 90
    ),
    plot.subtitle = element_text(
      margin = margin(10,0,0,0, "pt"),
      hjust = 0.5,
      size = 70,
      lineheight = 0.3
    ),
    plot.caption = element_textbox(
      hjust = 0.5,
      halign = 0.5
    ),
    panel.grid = element_line(
      colour = "grey80",
      linewidth = 0.1
    ),
    legend.title.position = "top",
    legend.margin = margin(0,0,0,0, "pt"),
    legend.box.margin = margin(0,0,0,0, "pt"),
    legend.text = element_text(
      margin = margin(0,0,0,2, "pt"),
      size = 60
    ),
    legend.title = element_text(
      margin = margin(0,0,5,0, "pt"),
      size = 60
    ),
    legend.background = element_rect(
      fill = "transparent",
      colour = "transparent"
    ),
    legend.key.height = unit(40, "pt")
  )

ggsave(
  plot = g,
  filename = here::here(
    "geocomputation", "images",
    "osm_packages_4.png"
  ),
  height = 3800,
  width = 2800,
  units = "px",
  bg = "white"
)```

Nice example, thanks! Would it be possible to make it more reproducible sharing how you created df1, wid1 and wid2? Thanks