#Setting up the environment
rm(list = ls()) #cleaning up the code environment
graphics.off() # cleans all the plots
#Loading required libraries, downloading required data from online sources, and loading the data:
library(sf) # simple features' library
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(spData) # library of spatial datasets
## 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, ...
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ 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
library(ggplot2)
library(readxl)
#Ethipia admin boundaries (Level 3- woredas): https://data.humdata.org/dataset/cod-ab-eth
#population density by woredas: https://data.humdata.org/dataset/cod-ps-eth
#major roads: https://data.humdata.org/dataset/hotosm_eth_roads
#high voltage electricity grid: https://datacatalog.worldbank.org/search/dataset/0039311
#ERSS enumeration areas (rural; surveyed twice): https://microdata.worldbank.org/index.php/catalog/2053
setwd("/Users/vaibhavjain/Desktop/BSE/Term 2/1. Geospatial/Assignment 1/Paper 2_Ethiopia")
# Load the world map
gmsf.world <- world
# Filter for Ethiopia's geometry
sf.Ethiopia <- gmsf.world%>% filter(iso_a2=="ET")
#reading woreda shapefiles
sf.Ethiopia_subadmin <- st_read("eth_admbnda_adm3_csa_bofedb_2021.shp")
## Reading layer `eth_admbnda_adm3_csa_bofedb_2021' from data source
## `/Users/vaibhavjain/Desktop/BSE/Term 2/1. Geospatial/Assignment 1/Paper 2_Ethiopia/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
#map with admin boundaries
ggplot(data = sf.Ethiopia) +
geom_sf(data = sf.Ethiopia_subadmin, color = "grey") + # Ethiopian woredas
coord_sf(crs = st_crs(4326), expand = FALSE, datum = NA) + # Force flat projection
theme_minimal() +
labs(title = "Ethiopia: Woredas (Level 3 Subnational Boundaries)") +
theme(
panel.grid = element_blank() # Remove gridlines entirely
)
#opening the population file
pop_data <- read_excel("eth_admpop_2023.xlsx", sheet = 5)
#renaming ID column for merge & executing the left-join merge
pop_data <- pop_data %>%
rename(ADM3_PCODE = admin3Pcode)
sf.Ethiopia_subadmin <- sf.Ethiopia_subadmin %>%
left_join(pop_data %>% select(ADM3_PCODE, Total), by = "ADM3_PCODE")
#population density
sf.Ethiopia_subadmin <- sf.Ethiopia_subadmin %>%
mutate(population_density = Total / sum(Total, na.rm = TRUE))
#map with population density and admin boundaries
ggplot(data = sf.Ethiopia) +
geom_sf(data = sf.Ethiopia_subadmin, aes(fill = population_density), color = "grey") + # Ethiopian woredas
scale_fill_gradient(name = "Population Density", low = "#d0e7f2", high = "#2a5674") +
coord_sf(crs = st_crs(4326), expand = FALSE, datum = NA) + # Force flat projection
theme_minimal() +
labs(title = "Ethiopia: Population Density by Woredas (Level 3 Subnational Boundaries)") +
theme(
panel.grid = element_blank() # Remove gridlines entirely
)
#reading major roads shapefiles
sf.Ethiopia_roads <- st_read("hotosm_eth_roads_lines_shp/hotosm_eth_roads_lines_shp.shp")
## Reading layer `hotosm_eth_roads_lines_shp' from data source
## `/Users/vaibhavjain/Desktop/BSE/Term 2/1. Geospatial/Assignment 1/Paper 2_Ethiopia/hotosm_eth_roads_lines_shp/hotosm_eth_roads_lines_shp.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 420912 features and 18 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 32.9902 ymin: 3.410445 xmax: 47.98158 ymax: 14.8554
## Geodetic CRS: WGS 84
unique(sf.Ethiopia_roads$highway)
## [1] "tertiary" "residential" "unclassified" "motorway_link"
## [5] "primary" "service" "track" "escape"
## [9] "path" "secondary" "footway" "construction"
## [13] "trunk_link" "steps" "primary_link" "trunk"
## [17] "secondary_link" "tertiary_link" "road" "living_street"
## [21] "motorway" "bridleway" "pedestrian" "platform"
## [25] "corridor" "services" "proposed" "emergency_bay"
# Define major highways
major_highways <- c("motorway", "trunk")
# Filter the roads dataset
sf.Ethiopia_roads <- sf.Ethiopia_roads %>%
filter(highway %in% major_highways)
unique(sf.Ethiopia_roads$highway)
## [1] "trunk" "motorway"
#map with population density, admin boundaries, and major roads
ggplot(data = sf.Ethiopia) +
geom_sf(data = sf.Ethiopia_subadmin, aes(fill = population_density), color = "grey") + # Ethiopian woredas
geom_sf(data = sf.Ethiopia_roads, color = "black", size = 0.3) + # Add major roads
scale_fill_gradient(name = "Population Density", low = "#d0e7f2", high = "#2a5674") +
coord_sf(crs = st_crs(4326), expand = FALSE, datum = NA) + # Force flat projection
theme_minimal() +
labs(title = "Ethiopia: Population Density by Woredas (Level 3 Subnational Boundaries) & Major Roads") +
theme(
panel.grid = element_blank() # Remove gridlines entirely
)
#reading electricity transmission shapefiles
sf.Ethiopia_elec <- st_read("ethiopia-electricity-transmission-network/Ethiopia Electricity Transmission Network.shp")
## Reading layer `Ethiopia Electricity Transmission Network' from data source
## `/Users/vaibhavjain/Desktop/BSE/Term 2/1. Geospatial/Assignment 1/Paper 2_Ethiopia/ethiopia-electricity-transmission-network/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
unique(sf.Ethiopia_elec$VOLTAGE_KV)
## [1] 66 132 230 330 45 400
# Define major transmission lines
major_lines <- c("230", "330", "400")
# Filter the roads dataset
sf.Ethiopia_elec <- sf.Ethiopia_elec %>%
filter(VOLTAGE_KV %in% major_lines)
unique(sf.Ethiopia_elec$VOLTAGE_KV)
## [1] 230 330 400
#map with population density, admin boundaries, major roads, and major electricity transmission lines
ggplot(data = sf.Ethiopia) +
geom_sf(data = sf.Ethiopia_subadmin, aes(fill = population_density), color = "grey") + # Ethiopian woredas
geom_sf(data = sf.Ethiopia_roads, color = "black", size = 0.3) + # Add major roads
geom_sf(data = sf.Ethiopia_elec, color = "red", size = 0.5) + # Add electricity infrastructure in red
scale_fill_gradient(name = "Population Density", low = "#d0e7f2", high = "#2a5674") +
coord_sf(crs = st_crs(4326), expand = FALSE, datum = NA) + # Force flat projection
theme_minimal() +
labs(title = "Ethiopia: Population Density by Woredas (Level 3 Subnational Boundaries), Major Roads, and Major Transmission Lines") +
theme(
panel.grid = element_blank() # Remove gridlines entirely
)
#ERSS Survey coordinates
ERSS <- read.csv("pub_eth_householdgeovariables_y1.csv")
ERSS_sf <- st_as_sf(ERSS, coords = c("LON_DD_MOD", "LAT_DD_MOD"), crs = 4326)
#map with population density, admin boundaries, major roads, and major electricity transmission lines
ggplot(data = sf.Ethiopia) +
geom_sf(data = sf.Ethiopia_subadmin, aes(fill = population_density), color = "grey") +
geom_sf(data = sf.Ethiopia_roads, aes(color = "Major Roads"), size = 0.3) +
geom_sf(data = sf.Ethiopia_elec, aes(color = "Transmission Lines"), size = 0.5) +
geom_sf(data = ERSS_sf, aes(color = "ERSS Locations"), size = 1, alpha = 0.6) +
# Adjust fill gradient for population density
scale_fill_gradient(name = "Population Density", low = "#d0e7f2", high = "#2a5674") +
# Add a custom legend for colors
scale_color_manual(
name = "Legend",
values = c(
"Major Roads" = "black", # Black for roads
"Transmission Lines" = "red", # Red for transmission lines
"ERSS Locations" = "black" # Blue for points
)
) +
coord_sf(crs = st_crs(4326), expand = FALSE, datum = NA) +
labs(title = "Ethiopian Population Density and ERSS Sample Villages with Major Roads & Electricity Transmission Lines") +
theme_minimal() +
theme(
panel.grid = element_blank() # Remove gridlines entirely
)
