#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(rnaturalearth)
library(ggplot2)
library(geosphere)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(spData)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(readxl)
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(viridis)
## Loading required package: viridisLite
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
#load the world dataset
sf.world <- world
sf.world2<- ne_countries()
# Add the population estimate (pop_est) from the NE dataset to the world dataset
sf.world <- sf.world %>%
left_join(st_drop_geometry(sf.world2)[, c("iso_a2", "pop_est")], by = "iso_a2")
# Renaming the columns such that 'pop_est' is 'pop'
sf.world$pop_temp <- sf.world$pop
sf.world$pop <- NULL
sf.world$pop <- sf.world$pop_est
sf.world$pop_est <- NULL
# Compute centroids for sf.world and extract coordinates
sf.world.centroids <- st_centroid(sf.world) %>%
st_coordinates() %>%
as.data.frame()
## Warning: st_centroid assumes attributes are constant over geometries
# Merge the coordinates with population data
sf.world.pop <- cbind(sf.world, sf.world.centroids)
# Modify sf.world to include a log-transformed population column
sf.world$pop_log <- log10(sf.world$pop)
# Plot
ggplot(data = sf.world) +
# Use fill color to represent population density (log scale)
geom_sf(aes(fill = pop_log), color = "black", size = 0.3) + # Black outlines for borders
# Use a blue color palette
scale_fill_viridis_c(option = "plasma", name = "Log Population", direction = -1) +
# Add title and labels
labs(
title = "World Population Distribution",
subtitle = "Countries shaded by population (log scale)",
caption = "Source: Natural Earth Data"
) +
# Improve theme aesthetics
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 16), # Make title bold
plot.subtitle = element_text(size = 12),
legend.position = "bottom",
panel.background = element_rect(fill = "white", color = NA),
axis.title = element_blank(),
axis.text = element_blank(),
panel.grid.major = element_blank()
)
# Read the airports shapefile
sf.airports <- st_read("/Users/vaibhavjain/Downloads/ne_50m_airports/ne_50m_airports.shp")
## Reading layer `ne_50m_airports' from data source
## `/Users/vaibhavjain/Downloads/ne_50m_airports/ne_50m_airports.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 281 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -175.1356 ymin: -53.00507 xmax: 178.56 ymax: 71.2893
## Geodetic CRS: WGS 84
sf.airports <- st_transform(sf.airports, 4326)
# Filter for terminals
sf.airports <- sf.airports %>% filter(location == "terminal")
# Extract coordinates for sf.airports
sf.airports.coords <- st_coordinates(sf.airports) %>%
as.data.frame()
# Merge coordinates back into sf.airports for plotting
sf.airports <- cbind(sf.airports, sf.airports.coords)
# Plot
ggplot(data = sf.world) +
# Use fill color to represent population density (log scale)
geom_sf(aes(fill = pop_log), color = "black", size = 0.3) + # Black outlines for borders
# Use a blue color palette
scale_fill_viridis_c(option = "plasma", name = "Log Population", direction = -1) +
# Plot airports as blue points
geom_point(data = sf.airports, aes(x = X, y = Y), color = "yellow", size = 0.25, alpha = 0.8) +
# Add title and labels
labs(
title = "World Population Distribution with Airports",
subtitle = "Countries shaded by population (log scale), airports marked in yellow",
caption = "Source: Natural Earth Data"
) +
# Improve theme aesthetics
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 16), # Make title bold
plot.subtitle = element_text(size = 12),
legend.position = "bottom",
panel.background = element_rect(fill = "white", color = NA),
axis.title = element_blank(),
axis.text = element_blank(),
panel.grid.major = element_blank()
)
# Read the ports shapefile
sf.ports <- st_read("/Users/vaibhavjain/Downloads/ne_10m_ports/ne_10m_ports.shp")
## Reading layer `ne_10m_ports' from data source
## `/Users/vaibhavjain/Downloads/ne_10m_ports/ne_10m_ports.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1081 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -171.758 ymin: -54.80944 xmax: 179.3094 ymax: 78.22611
## Geodetic CRS: WGS 84
sf.ports <- st_transform(sf.ports, 4326)
# Extract coordinates for sf.ports
sf.ports.coords <- st_coordinates(sf.ports) %>%
as.data.frame()
# Merge coordinates back into sf.airports for plotting
sf.ports <- cbind(sf.ports, sf.ports.coords)
# Plot
ggplot(data = sf.world) +
# Use fill color to represent population density (log scale)
geom_sf(aes(fill = pop_log), color = "black", size = 0.3) + # Black outlines for borders
# Use a blue color palette
scale_fill_viridis_c(option = "plasma", name = "Log Population", direction = -1) +
# Plot airports as blue points
geom_point(data = sf.airports, aes(x = X, y = Y), color = "yellow", size = 0.25, alpha = 0.8) +
# Plot ports as orange points
geom_point(data = sf.ports, aes(x = X, y = Y), color = "green", size = 0.25, alpha = 0.8) +
# Add title and labels
labs(
title = "World Population Distribution with Airports and Ports",
subtitle = "Countries shaded by population (log scale), airports in yellow, ports in green",
caption = "Source: Natural Earth Data"
) +
# Improve theme aesthetics
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 16), # Make title bold
plot.subtitle = element_text(size = 12),
legend.position = "bottom",
panel.background = element_rect(fill = "white", color = NA),
axis.title = element_blank(),
axis.text = element_blank(),
panel.grid.major = element_blank()
)
# Histogram: Country Population by Continent
ggplot(sf.world, aes(x = pop)) +
# Improved histogram with color gradient
geom_histogram(aes(fill = ..count..), bins = 30, color = "black", alpha = 0.8) +
# Facet by continent
facet_wrap(~continent, scales = "free_y") +
# Log-scale for better distribution
scale_x_log10() +
# Use a color gradient for fill
scale_fill_viridis_c(option = "plasma", name = "Frequency") +
# Improved labels
labs(
title = "Country Population Distribution by Continent",
subtitle = "Histogram of country populations (log scale)",
x = "Population (log scale)",
y = "Frequency"
) +
# Improve theme aesthetics
theme_light() +
theme(
plot.title = element_text(face = "bold", size = 16, color = "black"),
plot.subtitle = element_text(size = 12, color = "gray30"),
strip.text = element_text(face = "bold", size = 12), # Bold facet titles
legend.position = "right",
panel.grid.minor = element_blank() # Remove minor grid lines for cleaner look
)
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).
# Download populated places dataset
populated_places <- ne_download(scale = "medium", type = "populated_places",
category = "cultural", returnclass = "sf")
## Reading layer `ne_50m_populated_places' from data source
## `/private/var/folders/yx/_q0vqmws1396l1rbtrq0nn7h0000gn/T/RtmpoGlvWK/ne_50m_populated_places.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1251 features and 137 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -175.2206 ymin: -90 xmax: 179.2166 ymax: 78.22097
## Geodetic CRS: WGS 84
#Ensure CRS consistency
populated_places <- st_transform(populated_places, 4326)
# Select top 20 most populated cities per country
top_cities <- populated_places %>%
arrange(desc(POP_MAX)) %>% # Sort by largest population first
group_by(ADM0_A3) %>% # Group by country (ISO3 code)
slice_head(n = 20) %>% # Select top 20 cities per country
ungroup()
#finding out countries that are too small/overseas dependent territories etc. and are not there in our main dataset (sf.world)
extra_values <- setdiff(unique(top_cities$ISO_A2), unique(sf.world$iso_a2))
#we remove the extra_values from top_cities so that there aren't errors later while calculating distances
top_cities <- top_cities %>%
filter(!ISO_A2 %in% extra_values)
#For Euclidean distance, we would have to change the data to a different projection. However, since we already have the latitude and longitude of the data available, we can calculate the more accurate Great-Circle Distance
# Merge ports and airports into a single dataset for unified distance calculations
key_locations <- bind_rows(sf.ports, sf.airports)
# Perform a spatial join to add continent info to top_cities
top_cities <- st_join(top_cities, sf.world %>% select(continent, iso_a2), join = st_within)
# Compute average distance to the nearest transport hub (port or airport) per country
distance_data <- top_cities %>%
filter(!st_is_empty(geometry)) %>%
group_by(ADM0_A3, ADM0NAME, continent) %>% # Group by country
summarise(
avg_distance = mean(sapply(1:n(), function(i) {
if (is.na(geometry[i])) {
print(paste("DEBUG: Skipping empty geometry for", ADM0_A3[i], ADM0NAME[i]))
return(NA)
}
min_dist <- min(st_distance(st_geometry(top_cities)[i], st_geometry(key_locations), by_element = FALSE), na.rm = TRUE)
return(min_dist)
})), .groups = "drop"
)
# Convert meters to kilometers
distance_data <- distance_data %>%
mutate(avg_distance = avg_distance / 1000)
#removing NA data
distance_data <- distance_data %>% filter(!is.na(continent))
#Histogram for Average Distance per Country by Continent
ggplot(distance_data, aes(x = avg_distance)) +
# Improved histogram with color gradient
geom_histogram(aes(fill = ..count..), bins = 30, color = "black", alpha = 0.8) +
# Facet by continent
facet_wrap(~continent, scales = "free_y") +
# Log-scale for better distribution
scale_x_log10() +
# Use a color gradient for fill
scale_fill_viridis_c(option = "plasma", name = "Frequency") +
# Improved labels
labs(
title = "Average Distance from Cities to Nearest Port/Airport by Continent",
subtitle = "Histogram of country-level average distances (log scale)",
x = "Average Distance (km, log scale)",
y = "Frequency"
) +
# Improve theme aesthetics
theme_light() +
theme(
plot.title = element_text(face = "bold", size = 16, color = "black"),
plot.subtitle = element_text(size = 12, color = "gray30"),
strip.text = element_text(face = "bold", size = 12), # Bold facet titles
legend.position = "right",
panel.grid.minor = element_blank() # Remove minor grid lines for cleaner look
)
#Loading the relevant shapefile
Market_locations <- st_read("/Users/vaibhavjain/Downloads/116359-V1/data/1.-Price--Production--and-Population-Data/Catchments/MktCatch6Plus5_Dissolve.shp")
## Reading layer `MktCatch6Plus5_Dissolve' from data source
## `/Users/vaibhavjain/Downloads/116359-V1/data/1.-Price--Production--and-Population-Data/Catchments/MktCatch6Plus5_Dissolve.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 235 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -1690159 ymin: -4179878 xmax: 4948778 ymax: 3356116
## Projected CRS: World_Behrmann
plot(Market_locations)
#loading the market location data
market_centroids <- read_excel("/Users/vaibhavjain/Downloads/116359-V1/data/1.-Price--Production--and-Population-Data/MktCoords.xlsx")
#computing market centroids
market_centroids_sf <- market_centroids %>%
st_as_sf(coords = c("longitude", "latitude"), crs=4326)
#africa base data
africa <- ne_countries(continent = "Africa", scale = "medium", returnclass = "sf")
# Creating the Plot
ggplot() +
# Base Map with Elegant Styling
geom_sf(data = africa, fill = "#f0f0f0", color = "#595959", linewidth = 0.4) +
# Market Centroids with Enhanced Visibility
geom_sf(data = market_centroids_sf,
color = "#d73027", size = 2, alpha = 0.85) +
# Custom Titles and Labels
labs(title = "Market Locations Across Africa",
subtitle = "Geospatial visualization of market centroids") +
# Elegant Theme Adjustments
theme_void() + # Removes grid and unnecessary elements
theme(
plot.background = element_rect(fill = "#ffffff", color = NA), # White background
panel.grid = element_blank(), # Remove grid
axis.text = element_blank(), # Remove axis labels
axis.ticks = element_blank(),
plot.title = element_text(face = "bold", size = 18, color = "#333333", hjust = 0.5),
plot.subtitle = element_text(size = 14, color = "#666666", hjust = 0.5),
plot.caption = element_text(size = 10, color = "#888888", hjust = 0.5)
)
#filtering africa from world
African_countries <- world%>%
filter(world$continent=="Africa")
#airports in Africa
airports_in_africa <- st_filter(sf.airports, African_countries)
#Plot the map
ggplot() +
# Base Map with Soft, Elegant Colors
geom_sf(data = African_countries, aes(fill = name_long),
color = "#595959", size = 0.3, alpha = 0.6, show.legend = FALSE) +
# Airports Layer with Better Visibility
geom_sf(data = airports_in_africa,
color = "#0047AB", size = 1.5, alpha = 0.8) +
# Custom Titles and Labels
labs(title = "Airports Located in Africa",
subtitle = "Filtered using spatial intersection") +
# Elegant Theme Adjustments
theme_void() + # Removes grid and unnecessary elements
theme(
plot.background = element_rect(fill = "#ffffff", color = NA), # Clean white background
panel.grid = element_blank(), # Remove grid
axis.text = element_blank(), # Remove axis labels
axis.ticks = element_blank(),
plot.title = element_text(face = "bold", size = 18, color = "#333333", hjust = 0.5),
plot.subtitle = element_text(size = 14, color = "#666666", hjust = 0.5),
plot.caption = element_text(size = 10, color = "#888888", hjust = 0.5)
)
#loading the roads data from natural earth
roads <- ne_download(
scale = 10,
type = "roads",
category = "cultural",
returnclass = "sf"
)
## Reading layer `ne_10m_roads' from data source
## `/private/var/folders/yx/_q0vqmws1396l1rbtrq0nn7h0000gn/T/RtmpoGlvWK/ne_10m_roads.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 56600 features and 31 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -166.5325 ymin: -55.11212 xmax: 178.4191 ymax: 71.17768
## Geodetic CRS: WGS 84
#loading the coastline data from natural earth
coastline <- ne_download(
scale = 10,
type = "coastline",
category = "physical",
returnclass = "sf"
)
## Reading layer `ne_10m_coastline' from data source
## `/private/var/folders/yx/_q0vqmws1396l1rbtrq0nn7h0000gn/T/RtmpoGlvWK/ne_10m_coastline.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 4133 features and 3 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -180 ymin: -85.22194 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84
#filtering the roads and coastline for just Africa
roads_in_africa <- st_filter(roads, African_countries)
coastline_africa <- st_filter(coastline, African_countries)
# Create the Map
ggplot() +
# Base Map with Elegant Styling
geom_sf(data = African_countries, aes(fill = name_long),
color = "#4A4A4A", size = 0.3, alpha = 0.5, show.legend = FALSE) +
# Roads with Subtle Contrast
geom_sf(data = roads_in_africa, color = "#9E9E9E", size = 0.6, alpha = 0.8) +
# Coastline for Emphasis
geom_sf(data = coastline_africa, color = "#00BFFF", size = 0.7, alpha = 0.9) +
# Airports Layer in Distinctive Deep Blue
geom_sf(data = airports_in_africa, color = "#0047AB", size = 2.5, alpha = 0.9) +
# Market Centroids in a Bold, Strong Red
geom_sf(data = market_centroids_sf, color = "#B22222", size = 1.5, alpha = 0.7) +
# Title, Labels & Caption Enhancements
labs(title = "Markets, Airports, Roads & Coastlines in Africa",
subtitle = "A spatial visualization of infrastructure & economic hubs") +
# Elegant Theme Adjustments
theme_void() + # Removes unnecessary elements for a clean look
theme(
plot.background = element_rect(fill = "#ffffff", color = NA), # White background
panel.grid = element_blank(), # Remove grid
axis.text = element_blank(), # Remove axis labels
axis.ticks = element_blank(),
plot.title = element_text(face = "bold", size = 20, color = "#222222", hjust = 0.5),
plot.subtitle = element_text(size = 14, color = "#555555", hjust = 0.5),
plot.caption = element_text(size = 10, color = "#777777", hjust = 0.5)
)
#calculating the distance between the market and the coast, nearest road, and nearest airpot,
market_centroids_sf$distance_coast<- apply(st_distance(market_centroids_sf, coastline_africa), 1, min)
market_centroids_sf$distance_road<- apply(st_distance(market_centroids_sf, roads_in_africa), 1, min)
market_centroids_sf$distance_airport<- apply(st_distance(market_centroids_sf, airports_in_africa), 1, min)
#loading the market prices dataset
market_prices <- read_excel("/Users/vaibhavjain/Downloads/116359-V1/data/1.-Price--Production--and-Population-Data/PriceMaster4GAMS.xlsx")
# Convert 0s to NA
market_prices[market_prices == 0] <- NA
# Convert 0s to NA in numeric columns only
numeric_cols <- names(market_prices)[sapply(market_prices, is.numeric)]
market_prices[numeric_cols] <- lapply(market_prices[numeric_cols], function(x) ifelse(x == 0, NA, x))
#Delete Gulf and Bangkok that are not present in geospatial data
market_prices <- market_prices[market_prices$market != "Bangkok", ]
market_prices <- market_prices[market_prices$market != "Gulf", ]
unique(market_centroids$market)
## [1] "Luanda" "Cotonou" "Malanville"
## [4] "Natitingou" "Parakou" "Gaborone"
## [7] "Bobo Dioulasso" "Dedougou" "Fada Ngourma"
## [10] "Ouagadougou" "Bujumbura" "Gitega"
## [13] "Muyinga" "Bamenda" "Douala"
## [16] "Garoua" "Yaounde" "Bambari"
## [19] "Bangassou" "Bangui" "Abeche"
## [22] "Moundou" "Ndjamena" "Sarh"
## [25] "Brazzaville" "Impfondo" "Pointe Noire"
## [28] "Abengourou" "Abidjan" "Bouake"
## [31] "Daloa" "Man" "Odienne"
## [34] "Bandundu" "Bukavu" "Bunia"
## [37] "Butembo" "Gbadolite" "Goma"
## [40] "Isiro" "Kalemie" "Kamina"
## [43] "Kananga" "Kikwit" "Kindu"
## [46] "Kinshasa" "Kisangani" "Kolwezi"
## [49] "Lubumbashi" "Matadi" "Mbandaka"
## [52] "Mbuji Mayi" "Tshikapa" "Uvira"
## [55] "Zongo" "Djibouti" "Asmara"
## [58] "Massawa" "Addis Ababa" "Awasa"
## [61] "Bahir Dar" "Dessie" "Dire Dawa"
## [64] "Gambela" "Gode" "Gondar"
## [67] "Jijiga" "Jimma" "Mekele"
## [70] "Nekemte" "Yabelo" "Libreville"
## [73] "Banjul" "Accra" "Bolgatanga"
## [76] "Ho" "Kumasi" "Sekondi Takoradi"
## [79] "Tamale" "Wa" "Conakry"
## [82] "Kankan" "Labe" "Nzerekore"
## [85] "Bissau" "Eldoret" "Garissa"
## [88] "Kisumu" "Lodwar" "Mandera"
## [91] "Mombasa" "Moyale" "Nairobi"
## [94] "Nakuru" "Wajir" "Maseru"
## [97] "Gbarnga" "Monrovia" "Blantyre"
## [100] "Karonga" "Lilongwe" "Mangochi"
## [103] "Mzuzu" "Bamako" "Gao"
## [106] "Kayes" "Mopti" "Segou"
## [109] "Sikasso" "Adel Bagrou" "Nouakchott"
## [112] "Tintane" "Beira" "Chimoio"
## [115] "Cuamba" "Lichinga" "Maputo"
## [118] "Maxixe" "Milange" "Nacala"
## [121] "Nampula" "Pemba" "Quelimane"
## [124] "Tete" "Xai Xai" "Katima Mulilo"
## [127] "Oshakati" "Swakopmund" "Windhoek"
## [130] "Agadez" "Arlit" "Diffa"
## [133] "Maradi" "Niamey" "Tahoua"
## [136] "Zinder" "Abuja" "Akure"
## [139] "Benin City" "Calabar" "Enugu"
## [142] "Gombe" "Ibadan" "Ilorin"
## [145] "Jos" "Kaduna" "Kano"
## [148] "Katsina" "Lagos" "Lokoja"
## [151] "Maiduguri" "Makurdi" "Port Harcourt"
## [154] "Sokoto" "Yola" "Butare"
## [157] "Gisenyi" "Kigali" "Dakar"
## [160] "Kaolack" "Saint Louis" "Tambacounda"
## [163] "Touba" "Ziguinchor" "Bo"
## [166] "Freetown" "Kabala" "Baidoa"
## [169] "Beledweyne" "Bosaso" "Galkayo"
## [172] "Garoowe" "Hargeisa" "Kismayo"
## [175] "Mogadishu" "Bor" "Juba"
## [178] "Malakal" "Rumbek" "Wau"
## [181] "Ad Damazin" "Al Fashir" "Al Qadarif"
## [184] "El Geneina" "El Obeid" "Kadugli"
## [187] "Kassala" "Khartoum" "Kosti"
## [190] "Nyala" "Port Sudan" "Mbabane"
## [193] "Arusha" "Bukoba" "Dar es Salaam"
## [196] "Dodoma" "Iringa" "Kigoma"
## [199] "Mbeya" "Mtwara" "Musoma"
## [202] "Mwanza" "Singida" "Songea"
## [205] "Sumbawanga" "Tabora" "Tanga"
## [208] "Kara" "Lome" "Arua"
## [211] "Gulu" "Jinja" "Kampala"
## [214] "Lira" "Masindi" "Mbarara"
## [217] "Chipata" "Kabwe" "Kasama"
## [220] "Kitwe" "Livingstone" "Lusaka"
## [223] "Mongu" "Solwezi" "Bulawayo"
## [226] "Harare" "Hwange" "Masvingo"
## [229] "Mutare" "Johannesburg"
unique(market_prices$market)
## [1] "Luanda" "Cotonou" "Malanville" "Natitingou"
## [5] "Parakou" "Gaborone" "BoboDioulasso" "Dedougou"
## [9] "FadaNgourma" "Ouagadougou" "Bujumbura" "Gitega"
## [13] "Muyinga" "Bamenda" "Douala" "Garoua"
## [17] "Yaounde" "Bambari" "Bangassou" "Bangui"
## [21] "Abeche" "Moundou" "Ndjamena" "Sarh"
## [25] "Brazzaville" "Impfondo" "PointeNoire" "Abengourou"
## [29] "Abidjan" "Bouake" "Daloa" "Man"
## [33] "Odienne" "Bandundu" "Bukavu" "Bunia"
## [37] "Butembo" "Gbadolite" "Goma" "Isiro"
## [41] "Kalemie" "Kamina" "Kananga" "Kikwit"
## [45] "Kindu" "Kinshasa" "Kisangani" "Kolwezi"
## [49] "Lubumbashi" "Matadi" "Mbandaka" "MbujiMayi"
## [53] "Tshikapa" "Uvira" "Zongo" "Djibouti"
## [57] "Asmara" "Massawa" "AddisAbaba" "Awasa"
## [61] "BahirDar" "Dessie" "DireDawa" "Gambela"
## [65] "Gode" "Gondar" "Jijiga" "Jimma"
## [69] "Mekele" "Nekemte" "Yabelo" "Libreville"
## [73] "Banjul" "Accra" "Bolgatanga" "Ho"
## [77] "Kumasi" "SekondiTakoradi" "Tamale" "Wa"
## [81] "Conakry" "Kankan" "Labe" "Nzerekore"
## [85] "Bissau" "Eldoret" "Garissa" "Kisumu"
## [89] "Lodwar" "Mandera" "Mombasa" "Moyale"
## [93] "Nairobi" "Nakuru" "Wajir" "Maseru"
## [97] "Gbarnga" "Monrovia" "Blantyre" "Karonga"
## [101] "Lilongwe" "Mangochi" "Mzuzu" "Bamako"
## [105] "Gao" "Kayes" "Mopti" "Segou"
## [109] "Sikasso" "AdelBagrou" "Nouakchott" "Tintane"
## [113] "Beira" "Chimoio" "Cuamba" "Lichinga"
## [117] "Maputo" "Maxixe" "Milange" "Nacala"
## [121] "Nampula" "Pemba" "Quelimane" "Tete"
## [125] "XaiXai" "KatimaMulilo" "Oshakati" "Swakopmund"
## [129] "Windhoek" "Agadez" "Arlit" "Diffa"
## [133] "Maradi" "Niamey" "Tahoua" "Zinder"
## [137] "Abuja" "Akure" "BeninCity" "Calabar"
## [141] "Enugu" "Gombe" "Ibadan" "Ilorin"
## [145] "Jos" "Kaduna" "Kano" "Katsina"
## [149] "Lagos" "Lokoja" "Maiduguri" "Makurdi"
## [153] "PortHarcourt" "Sokoto" "Yola" "Butare"
## [157] "Gisenyi" "Kigali" "Dakar" "Kaolack"
## [161] "SaintLouis" "Tambacounda" "Touba" "Ziguinchor"
## [165] "Bo" "Freetown" "Kabala" "Baidoa"
## [169] "Beledweyne" "Bosaso" "Galkayo" "Garoowe"
## [173] "Hargeisa" "Kismayo" "Mogadishu" "Johannesburg"
## [177] "Bor" "Juba" "Malakal" "Rumbek"
## [181] "Wau" "AdDamazin" "AlFashir" "AlQadarif"
## [185] "ElGeneina" "ElObeid" "Kadugli" "Kassala"
## [189] "Khartoum" "Kosti" "Nyala" "PortSudan"
## [193] "Mbabane" "Arusha" "Bukoba" "DaresSalaam"
## [197] "Dodoma" "Iringa" "Kigoma" "Mbeya"
## [201] "Mtwara" "Musoma" "Mwanza" "Singida"
## [205] "Songea" "Sumbawanga" "Tabora" "Tanga"
## [209] "Kara" "Lome" "Arua" "Gulu"
## [213] "Jinja" "Kampala" "Lira" "Masindi"
## [217] "Mbarara" "Chipata" "Kabwe" "Kasama"
## [221] "Kitwe" "Livingstone" "Lusaka" "Mongu"
## [225] "Solwezi" "Bulawayo" "Harare" "Hwange"
## [229] "Masvingo" "Mutare"
# Compute the average price for each row (excluding non-numeric columns)
market_prices <- market_prices %>%
mutate(avg_price = rowMeans(select(., `1`:`144`), na.rm = TRUE)) # Adjust range as needed
market_prices <- market_prices[, -c(5:148)]
market_prices <- market_prices %>%
group_by(mktcode, country, market) %>%
summarise(avg_price = mean(avg_price, na.rm = TRUE), .groups = "drop")
market_prices <- left_join(market_centroids_sf, market_prices, by = c("mktcode", "market"))
# Create interactive scatter plots
p1 <- ggplot(market_prices, aes(x = distance_coast, y = avg_price, text = market)) +
geom_point(color = "#FF5733", size = 3, alpha = 0.8) +
labs(title = "Average Price vs Distance to Coast",
x = "Distance to Coast (km)",
y = "Average Price ($)") +
theme_minimal() +
theme(panel.background = element_rect(fill = "black"),
plot.background = element_rect(fill = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
text = element_text(color = "white"))
p2 <- ggplot(market_prices, aes(x = distance_road, y = avg_price, text = market)) +
geom_point(color = "#33FF57", size = 3, alpha = 0.8) +
labs(title = "Average Price vs Distance to Road",
x = "Distance to Road (km)",
y = "Average Price ($)") +
theme_minimal() +
theme(panel.background = element_rect(fill = "black"),
plot.background = element_rect(fill = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
text = element_text(color = "white"))
p3 <- ggplot(market_prices, aes(x = distance_airport, y = avg_price, text = market)) +
geom_point(color = "#337BFF", size = 3, alpha = 0.8) +
labs(title = "Average Price vs Distance to Airport",
x = "Distance to Airport (km)",
y = "Average Price ($)") +
theme_minimal() +
theme(panel.background = element_rect(fill = "black"),
plot.background = element_rect(fill = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
text = element_text(color = "white"))
# Convert ggplot to interactive plotly
p1_interactive <- ggplotly(p1, tooltip = "text")
p2_interactive <- ggplotly(p2, tooltip = "text")
p3_interactive <- ggplotly(p3, tooltip = "text")
# Print plots explicitly
p1_interactive
p2_interactive
p3_interactive