Geospatial

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

library(sf)

## Warning: package 'sf' was built under R version 4.4.2

## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE

library(ggplot2)
library(terra)

## Warning: package 'terra' was built under R version 4.4.2

## terra 1.8.15

# Set working directory
setwd("D:\\BSE Semester 2\\Geospatial Data Science and Economic Spatial Models\\Assignments\\Assign 1\\Question 1\\Q1")


# Question 1

# Description: To replicate the figure, we use data from different sources:
# 1. Global Dam Watch (GDW): to extract dam locations in South Africa. As in the paper, we filter the dataset to include only dams constructed up to 2002. (link: https://www.globaldamwatch.org/database) 
# 2. Global Administrative Areas (GADM): to extract district boundaries in South Africa. (link: https://gadm.org/download_country.html)
# 3. RCMRD Africa GeoPortal: as a proxy for river gradients, we calculate the average land slope for each district using elevation data (link: https://rcmrd.africageoportal.com/datasets/aafd352ca9044f23a78274e782dfe788/explore)


# Set working directory
setwd("D:\\BSE Semester 2\\Geospatial Data Science and Economic Spatial Models\\Assignments\\Assign 1\\Question 1\\Q1")

# Load Global Dam Watch (GDW) dataset containing dam locations
represas <- st_read("D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Question 1/Q1/GDW_reservoirs_v1_0.shp")

## Reading layer `GDW_reservoirs_v1_0' from data source 
##   `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Question 1\Q1\GDW_reservoirs_v1_0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 35295 features and 71 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -159.4638 ymin: -45.88067 xmax: 178.1558 ymax: 70.39624
## Geodetic CRS:  WGS 84

# Filter dams located in South Africa and constructed before the end of 2002
represas_sa <- represas[represas$COUNTRY == "South Africa" & represas$YEAR_DAM <= 2002, ]

# Load GADM dataset containing district boundaries for South Africa
distritos <- st_read("gadm41_ZAF_3.shp")

## Reading layer `gadm41_ZAF_3' from data source 
##   `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Question 1\Q1\gadm41_ZAF_3.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 234 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 16.45189 ymin: -34.83514 xmax: 32.89125 ymax: -22.12503
## Geodetic CRS:  WGS 84

# Load Digital Elevation Model (DEM) - Terrain elevation data
srtm <- rast("South_Africa_SRMT30meters.tif")

# Reduce the resolution of the DEM raster 
srtm_reducido <- aggregate(srtm, fact = 5, fun = mean)

## |---------|---------|---------|---------|=========================================

# Calculate terrain slope in degrees from the DEM
slope <- terrain(srtm_reducido, v = "slope", unit = "degrees")

# Extract the average slope value for each district
gradiente_por_distrito <- extract(slope, vect(distritos), fun = mean, na.rm = TRUE)

# Assign the calculated average slope to the district shapefile
distritos$avg_slope <- gradiente_por_distrito[, 2]

# plot dam locations and terrain slope by district
ggplot() +
  # districts, colored by average slope
  geom_sf(data = distritos, aes(fill = avg_slope), color = "black") +
  # dam locations
  geom_sf(data = represas_sa, aes(color = "Dam location"), size = 1) +
  # Inverted color scale for slope
  scale_fill_viridis_c(option = "plasma", name = "Average Slope", direction = -1) +
  # Color scale for dams
  scale_color_manual(values = c("Dam location" = "blue"), name = "Legend") +
  # Labels and title
  labs(title = "Dam Locations and Terrain Slope by District",
       caption = "Source: GADM, RCMRD & GDW")

# Question 2

# -------------------------------------------------------------------
# Title: Population Density, Roads, and Electricity in Ethiopia
# Data Source:
# -------------------------------------------------------------------
# Data Sources:
# 1. Administrative Boundaries (Ethiopia, ADM Level 3)  
#    Source: Ethiopian Central Statistical Agency (CSA), 2021  
#    File: eth_admbnda_adm3_csa_bofedb_2021.shp  

# 2. Population Density Data  
#    Source: Humanitarian Data Exchange (HumData), 2022  
#    URL: https://data.humdata.org/dataset/cod-ps-eth  
#    File: pop_adm3_2022_v2.csv  

# 3. Ethiopia Road Network Data  
#    Source: Ethiopian Road Authority (ERA), 2022  
#    File: Ethiopia_Roads.shp  

# 4. Ethiopia Electricity Transmission Grid  
#    Source: Ethiopian Electric Power Corporation (EEP), 2022  
#    File: Ethiopia Electricity Transmission Network.shp  

# 5. Enumeration Localities (Survey Data)  
#    Source: World Bank - Ethiopian Rural Socioeconomic Survey (ERSS), 2022  
#    URL: https://microdata.worldbank.org/index.php/catalog/2053/data-dictionary  
#    File: pub_eth_householdgeovariables_y1.csv  

# 6. Ethiopia Power Plants Data  
#    Source: World Bank - Ethiopia Power Plants Dataset, 2022  
#    URL: https://datacatalog.worldbank.org/search/dataset/0041714/Ethiopia-Power-Plants  
#    File: ETH_PowerPlants.shp  
# -------------------------------------------------------------------
# Paper: Fried, B., Lagakos, D., Poschke, M., & Santangelo, G. (2021). 
#        "The role of electricity in structural transformation: Evidence from Ethiopia."
#        *Journal of Urban Economics.* https://doi.org/10.1016/j.jue.2020.103293
# -------------------------------------------------------------------

# install.packages("ggrepel")

# Load required libraries
library(sf)            # For handling spatial data
library(dplyr)         # For data manipulation

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:terra':
## 
##     intersect, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)       # For visualization
library(ggrepel)       # For improved text labels

## Warning: package 'ggrepel' was built under R version 4.4.2

library(tidyverse)     # For data transformation

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ✔ readr     2.1.5

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# -------------------------------------------------------------------
# Set working directory and define paths for GIS data
# -------------------------------------------------------------------
setwd("D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Questions 2 and 5/Assignment_1/q2/p2")
path <- getwd()



# -------------------------------------------------------------------
# Data Loading and Preprocessing
# -------------------------------------------------------------------

# 1. Load Administrative Boundaries (Ethiopia, ADM Level 3)
sf.eth_admin_3 <- st_read(file.path(path, "adm_3/eth_admbnda_adm3_csa_bofedb_2021.shp"))

## Reading layer `eth_admbnda_adm3_csa_bofedb_2021' from data source 
##   `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Questions 2 and 5\Assignment_1\q2\p2\adm_3\eth_admbnda_adm3_csa_bofedb_2021.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1082 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 32.9918 ymin: 3.40667 xmax: 47.98824 ymax: 14.84548
## Geodetic CRS:  WGS 84

# 2. Load Population Density Data
eth_pop <- read.csv(file.path(path, "pop_adm3_2022_v2.csv"))

# Ensure column names match for merging
if (!("ADM3_PCODE" %in% colnames(eth_pop))) {
  eth_pop <- rename(eth_pop, ADM3_PCODE = admin3Pcode) 
}

# Join Population Data with Administrative Boundaries
joined_pop_geo <- left_join(sf.eth_admin_3, eth_pop, by = "ADM3_PCODE")

# 3. Load Main Roads Data
sf.roads <- st_read(file.path(path, "roads_2/Ethiopia_Roads.shp"))

## Reading layer `Ethiopia_Roads' from data source 
##   `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Questions 2 and 5\Assignment_1\q2\p2\roads_2\Ethiopia_Roads.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 76 features and 16 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 35.58148 ymin: 5.6777 xmax: 43.5421 ymax: 14.33759
## Geodetic CRS:  WGS 84

# Convert roads data to spatial format and remove certain links
roads_geo <- st_as_sf(sf.roads, wkt = "geometry", crs = 4326) %>% 
  filter(LINKNO != "A5-5")

# 4. Load Electricity Transmission Lines Data
sf.lines <- st_read(file.path(path, "electricity_grid/Ethiopia Electricity Transmission Network.shp"))

## Reading layer `Ethiopia Electricity Transmission Network' from data source 
##   `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Questions 2 and 5\Assignment_1\q2\p2\electricity_grid\Ethiopia Electricity Transmission Network.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 8 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 34.5333 ymin: 3.613251 xmax: 43.58116 ymax: 14.27694
## Geodetic CRS:  WGS 84

# Convert electricity grid data to spatial format
lines_geo <- st_as_sf(sf.lines, wkt = "geometry", crs = 4326)

# 5. Load Enumeration Localities Data (Survey Data)
enum <- read.csv(file.path(path, "erss_survey/pub_eth_householdgeovariables_y1.csv"))

# 6. Load Power Plants Data
plants <- st_read(file.path(path, "eth_powerplants/ETH_PowerPlants.shp"))

## Reading layer `ETH_PowerPlants' from data source 
##   `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Questions 2 and 5\Assignment_1\q2\p2\eth_powerplants\ETH_PowerPlants.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 22 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 36.9127 ymin: 6.856073 xmax: 40.17 ymax: 13.3
## CRS:           NA

# Convert Power Plants Data to Correct CRS
plants_geo <- st_set_crs(plants, 4326)

# -------------------------------------------------------------------
# Data Transformation
# -------------------------------------------------------------------

# Add a new column for Population Density (Log Scale)
joined_pop_geo$pupulation_density <- log(joined_pop_geo$T_TL / joined_pop_geo$Shape_Area)

# -------------------------------------------------------------------
# Plotting the Map of Ethiopia's Population Density and Infrastructure
# -------------------------------------------------------------------

ggplot() +
  # 1. Plot Population Density by Administrative Boundary
  geom_sf(data = joined_pop_geo, aes(fill = pupulation_density)) +
  scale_fill_gradient(low = "aliceblue", high = "midnightblue", name = "Population Density") +
  
  # 2. Overlay Main Roads (Dashed Black Lines)
  geom_sf(data = roads_geo, color = "black", size = 1, linetype = "dashed") +
  
  # 3. Overlay Electricity Transmission Lines (Red)
  geom_sf(data = lines_geo, color = "red", size = 0.5) +
  
  # 4. Plot Enumeration Localities (Small White Dots with Black Outline)
  geom_point(data = enum, aes(x = LON_DD_MOD, y = LAT_DD_MOD), 
             color = "black", fill = "white", size = 1, shape = 21) +
  
  # 5. Plot Power Plants (Blue Star Shape)
  geom_point(data = plants_geo, aes(x = LON, y = LAT), 
             shape = 5, color = "blue") +
  
  # 6. Add Title and Captions
  labs(
    title = "Population Density, Roads, and Electricity in Ethiopia",
    subtitle = "Overlay of Population Density with Roads, Electricity Lines, and Power Plants",
    caption = expression(bold("Figure 2: Map of Ethiopian Population Density, Infrastructure, and Power Grid"))
  ) +
  
  # 7. Style Adjustments for Clarity
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    plot.caption = element_text(size = 12, face = "bold", hjust = 0.5), # Centered Caption
    legend.position = "right"
  )

# Question 3

#Set working directory
# setwd("D:\\BSE Semester 2\\Geospatial Data Science and Economic Spatial Models\\Assignments\\Assign 1\\Question 3")


# Load necessary packages
library(sf)
library(geobr)

## Warning: package 'geobr' was built under R version 4.4.2

library(dplyr)
library(ggplot2)
library(readxl)
library(tidyr)
library(lwgeom)

## Warning: package 'lwgeom' was built under R version 4.4.2

## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.12.2, PROJ 9.4.1

## 
## Attaching package: 'lwgeom'

## The following object is masked from 'package:sf':
## 
##     st_perimeter

library(sfheaders)

## Warning: package 'sfheaders' was built under R version 4.4.2

# Read the .xls file with the Brazilian population for 2009
df <- read_excel("D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Question 3/UF_Municipio.xls", skip = 4, col_names = TRUE)

## New names:
## • `COD` -> `COD...2`
## • `COD` -> `COD...3`

# Get the geographical data from the package geobr
state <- read_state(year = 2010)

## Using year/date 2010

meso <- read_meso_region(year = 2010)

## Using year/date 2010

region <- read_region(year = 2010)

## Using year/date 2010

# Convert ESTIMADA (=population) into a numerical variable and handle missing values
# by setting them to 0
df <- df %>%
  mutate(ESTIMADA = as.numeric(ESTIMADA)) %>%
  mutate(ESTIMADA = replace_na(ESTIMADA, 0))

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `ESTIMADA = as.numeric(ESTIMADA)`.
## Caused by warning:
## ! NAs introduced by coercion

# Group by the country code
df_grouped <- df %>%
  group_by(COD...2) %>%
  summarise(
    Total_Estimada = sum(ESTIMADA, na.rm = TRUE),
    Count = n()
  )

# Remove the last row of df_merged - in df the last row was the total
# and it would have been considered another meso region
df_grouped <- df_grouped %>%
  slice(-n())

# Merge the geographical data with the population data
df_merged <- state %>%
  left_join(df_grouped, by = c("code_state" = "COD...2"))

# Create the column pop_share
df_merged <- df_merged %>%
  mutate(pop_share = Total_Estimada / sum(Total_Estimada) * 100)

# Aggregate regions as the author's did to get the correct boundaries
region <- region %>%
  mutate(
    name_region = case_when(
      name_region %in% c("Centro Oeste", "Norte") ~ 2,
      TRUE ~ 1
    )
  ) %>%
  dplyr::select(name_region, geom) %>%
  group_by(name_region) %>%
  summarise(geom = st_union(geom)) %>%
  sfheaders::sf_remove_holes()

# Final plot
ggplot(data = df_merged) +
  geom_sf(aes(fill = pop_share), color = "white") +
  scale_fill_gradientn(
    colors = c("lightblue", "blue", "darkblue"),
    name = "Pop. Share",
    limits = c(0, 25)
  ) +
  geom_sf(data = region, aes(color = factor(name_region)), fill = NA, linewidth = 1) +
  scale_color_manual(values = c("black", "red")) +
  guides(color = "none")

Including Plots

You can also embed plots, for example:

# Question 4

# Step 1  Import libraries
library(sf) # simple features' library
library(spData) # library of spatial datasets

## Warning: package 'spData' was built under R version 4.4.2

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(tidyverse) # dplyr, ggplot

# Set working directory
setwd("D:\\BSE Semester 2\\Geospatial Data Science and Economic Spatial Models\\Assignments\\Assign 1\\Question 4")


# Step 2 - Load shape files. Note have used different data source to create the map
# Vietnam map, data source - Link : https://gadm.org/download_country.html#google_vignette
sf.vietnam <- st_read('D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Question 4/gadm41_VNM_shp/gadm41_VNM_0.shp')

## Reading layer `gadm41_VNM_0' from data source 
##   `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Question 4\gadm41_VNM_shp\gadm41_VNM_0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 102.1446 ymin: 8.381355 xmax: 109.4692 ymax: 23.39269
## Geodetic CRS:  WGS 84

# Vietnam Roads, data source - Link : https://data.humdata.org/dataset/viet-nam-roads
# Note data is available for 2015
sf.roads <- st_read('D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Question 4/vnm_rdsl_2015_osm/vnm_rdsl_2015_OSM.shp')

## Reading layer `vnm_rdsl_2015_OSM' from data source 
##   `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Question 4\vnm_rdsl_2015_osm\vnm_rdsl_2015_OSM.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 118138 features and 5 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 102.1562 ymin: 8.596338 xmax: 109.4536 ymax: 23.37756
## Geodetic CRS:  WGS 84

# Geo map initial
# ggplot() + geom_sf(data = sf.vietnam) + geom_sf(data= sf.roads, aes(color=type))

# unique(sf.roads$type)

# Step 3 - Club and reduce road types
sf.roads <- sf.roads %>%
  mutate(
    road_group = case_when(
      type %in% c("motorway", "motorway_link", "") ~ "Freeway",
      type %in% c("primary", "secondary", "tertiary") ~ "Major Roads",
      type %in% c("residential", "service", "track") ~ "Minor Roads",
      TRUE ~ "Other"
    )
  )

# Define custom colors for each road group
road_colors <- c(
  "Highways" = "red",
  "Major Roads" = "blue",
  "Minor Roads" = "green",
  "Other" = "gray"
)

# Step 4 - Geo map with custom road groups and colors
ggplot() +
  geom_sf(data = sf.vietnam, fill = "lightyellow", color = "black") + # Vietnam map
  geom_sf(data = sf.roads, aes(color = road_group)) + # Roads with grouped types
  scale_color_manual(values = road_colors, name = "Road Types") + # Custom colors
  theme_minimal() +
  labs(title = "Road Map of Vietnam (2015)", subtitle = "Road Types Categorized by Function") +
  theme(legend.position = "bottom", legend.title = element_text(size = 12))

# Question 5

# -------------------------------------------------------------------
# Title: Map of Straight-Line Instrument and Radial Highways in Brazil
# Data Source: American Economic Association Replication Package
# Paper: Morten, M. and Oliveira, J. (2019). "Transportation, Migration, and Development: 
#           Evidence from Brazil." *American Economic Journal: Applied Economics*.
#           https://www.aeaweb.org/articles?id=10.1257/app.20180487
# -------------------------------------------------------------------

# Load required libraries
library(sf)            # For handling spatial data
library(spData)        # Spatial datasets
library(foreign)       # For working with data in foreign formats
library(ggrepel)       # For nicely repelled text labels in ggplot
library(tidyverse)     # For data manipulation and visualization

# -------------------------------------------------------------------
# Define custom themes to remove axes and adjust legends for a cleaner map
# -------------------------------------------------------------------
no_axis <- theme(
  axis.title = element_blank(),
  axis.text = element_blank(),
  axis.ticks = element_blank()
)

no_axis_legend <- theme(
  axis.title = element_blank(),
  axis.text = element_blank(),
  axis.ticks = element_blank(),
  legend.position = "bottom" # Place legend at the bottom
)

# -------------------------------------------------------------------
# Set working directory and define paths for GIS data
# -------------------------------------------------------------------
setwd("D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Questions 2 and 5/Assignment_1/q5")


path <- getwd()
path_data <- file.path(path, "GIS_data")

# Define a list for potential looping (if needed)
loop_list <- c("uf1940", "meso")

# -------------------------------------------------------------------
# Data Loading and Preprocessing
# -------------------------------------------------------------------

# Suppress unwanted output temporarily
sink(tempfile()) 

# 1. Load shapefile for the location of each state and simplify geometry
states <- st_read(file.path(path_data, "uf1940/uf1940_prj.shp"))
sf.states <- rmapshaper::ms_simplify(states, keep = 0.01, keep_shapes = TRUE)

# 2. Load shapefile for roads (year 2000) and simplify geometry
year <- 2000
file_name <- file.path(path_data, "roads", as.character(year), paste0("highways_", year, "_prj.shp"))
all_highways <- st_read(file_name)
sf.roads <- rmapshaper::ms_simplify(all_highways, keep = 0.01, keep_shapes = TRUE)

# 3. Load shapefile for Minimum Spanning Tree (MST) and simplify geometry
mst_pie <- st_read(file.path(path_data, "mst/mst_pie_prj.shp"))
sf.mst <- rmapshaper::ms_simplify(mst_pie, keep = 0.01, keep_shapes = TRUE)

# 4. Load shapefile for cities in Brazil and extract coordinates
capital_cities <- st_read(file.path(path_data, "cities/brazil_capital_cities_prj.shp"))
cities_xy <- cbind(capital_cities, st_coordinates(capital_cities))

# Restore output
sink()

# -------------------------------------------------------------------
# Plotting the Map of Brazil's Highway Network (2000)
# -------------------------------------------------------------------

ggplot() +
  # Plot states
  geom_sf(data = sf.states, fill = "white", color = "grey90") +
  
  # Plot minimum spanning tree
  geom_sf(data = sf.mst, size = 0.8, linetype = "11", aes(color = "Minimum spanning tree"), show.legend = "line") +
  
  # Plot non-radial highways (dashed)
  geom_sf(data = sf.roads %>% filter(dm_anlys_p == 1 & dm_radial == 0), size = 0.6, linetype = "dashed", 
          aes(color = "Non-radial highways (2000)"), show.legend = "line") +
  
  # Plot radial highways (solid)
  geom_sf(data = sf.roads %>% filter(dm_anlys_p == 1 & dm_radial == 1), size = 0.9, 
          aes(color = "Radial highways (2000)"), show.legend = "line") +
  
  # Add points for cities
  geom_point(data = cities_xy, aes(x = X, y = Y)) +
  
  # Add labels for cities
  geom_text_repel(data = cities_xy, aes(x = X, y = Y, label = CITY_NAME)) +
  
  # Set minimal theme and remove axes
  theme_minimal() +
  no_axis +
  
  # Add color legend
  labs(color = " ") +
  
  # Define manual colors for different elements
  scale_color_manual(
    values = c("#777676", "#868686", "#565555"),
    guide = guide_legend(override.aes = list(linetype = c("11", "dashed", "solid")))
  ) +
  
  # Add a properly formatted caption below the figure
  labs(
    caption = expression(bold("Figure 1: Map of straight-line instrument and radial highways"))
  ) +
  
  # Adjust theme to align the caption properly
  theme(
    plot.caption = element_text(hjust = 0.5, size = 12, face = "bold") # Center-aligned and bold caption
  )

Geospatial_PS1

Angad

2025-01-28

R Markdown

Including Plots