# Libraries
library(sf)
library(tidyverse)
library(spData)
library(readxl)
library(tidyr)
library(ggplot2)\(\color{darkblue}{\text{Part 1}}\)
1. Combination of the world with more data
Combination with:
- Population (point) data (do not use rasters!)
- Ports, airports, etc.
1.1. World dataset
Loading the world dataset:
## Simple feature collection with 6 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
## Geodetic CRS: WGS 84
## # A tibble: 6 × 11
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
## 2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
## 3 EH Western S… Africa Africa Northern… Inde… 9.63e4 NA NA
## 4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
## 5 US United St… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
## 6 KZ Kazakhstan Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
## # ℹ 2 more variables: gdpPercap <dbl>, geom <MULTIPOLYGON [°]>We keep only the variables that interest us for this task:
gmsf_world_sel <- gmsf_world %>% select(iso_a2, name_long, continent, type, pop, geom)
head(gmsf_world_sel)## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
## Geodetic CRS: WGS 84
## # A tibble: 6 × 6
## iso_a2 name_long continent type pop geom
## <chr> <chr> <chr> <chr> <dbl> <MULTIPOLYGON [°]>
## 1 FJ Fiji Oceania Soverei… 8.86e5 (((-180 -16.55522, -179.…
## 2 TZ Tanzania Africa Soverei… 5.22e7 (((33.90371 -0.95, 31.86…
## 3 EH Western Sahara Africa Indeter… NA (((-8.66559 27.65643, -8…
## 4 CA Canada North America Soverei… 3.55e7 (((-132.71 54.04001, -13…
## 5 US United States North America Country 3.19e8 (((-171.7317 63.78252, -…
## 6 KZ Kazakhstan Asia Soverei… 1.73e7 (((87.35997 49.21498, 86…Note: just using the world dataset to produce a map of total population by country is not enough, as there are countries with missing values (see below).
1.2. Population data by location
Source of the population data: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-populated-places/
gmsf_population <- st_read("/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/geo_assig2/map_1/populated_places/ne_10m_populated_places.shp")## Reading layer `ne_10m_populated_places' from data source
## `/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/geo_assig2/map_1/populated_places/ne_10m_populated_places.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7342 features and 137 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -179.59 ymin: -90 xmax: 179.3833 ymax: 82.48332
## Geodetic CRS: WGS 84Again, we only select those variables that may interest us:
gmsf_population_sel <- gmsf_population %>%
select(FEATURECLA, NAME, NAMEASCII, SOV0NAME, ADM0NAME, ISO_A2, POP_MAX, POP_MIN, geometry)
head(gmsf_population_sel)## Simple feature collection with 6 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58.304 ymin: -34.538 xmax: 0.7890036 ymax: 9.261
## Geodetic CRS: WGS 84
## FEATURECLA NAME NAMEASCII SOV0NAME
## 1 Admin-1 capital Colonia del Sacramento Colonia del Sacramento Uruguay
## 2 Admin-1 capital Trinidad Trinidad Uruguay
## 3 Admin-1 capital Fray Bentos Fray Bentos Uruguay
## 4 Admin-1 capital Canelones Canelones Uruguay
## 5 Admin-1 capital Florida Florida Uruguay
## 6 Admin-1 capital Bassar Bassar Togo
## ADM0NAME ISO_A2 POP_MAX POP_MIN geometry
## 1 Uruguay UY 21714 21714 POINT (-57.83612 -34.46979)
## 2 Uruguay UY 21093 21093 POINT (-56.901 -33.544)
## 3 Uruguay UY 23279 23279 POINT (-58.304 -33.139)
## 4 Uruguay UY 19698 19698 POINT (-56.284 -34.538)
## 5 Uruguay UY 32234 32234 POINT (-56.215 -34.099)
## 6 Togo TG 61845 61845 POINT (0.7890036 9.261)1.3. Population data by country
Source of the population by country data: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/
gmsf_population_by_countries <- st_read("/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/geo_assig2/map_1/countries/ne_10m_admin_0_countries.shp")## Reading layer `ne_10m_admin_0_countries' from data source
## `/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/geo_assig2/map_1/countries/ne_10m_admin_0_countries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 258 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84gmsf_population_by_countries_sel <- gmsf_population_by_countries %>%
select(SOVEREIGNT, ADMIN, ISO_A2, CONTINENT, POP_EST, POP_YEAR, geometry)
head(gmsf_population_by_countries_sel)## Simple feature collection with 6 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -109.4537 ymin: -55.9185 xmax: 140.9776 ymax: 7.35578
## Geodetic CRS: WGS 84
## SOVEREIGNT ADMIN ISO_A2 CONTINENT POP_EST POP_YEAR
## 1 Indonesia Indonesia ID Asia 270625568 2019
## 2 Malaysia Malaysia MY Asia 31949777 2019
## 3 Chile Chile CL South America 18952038 2019
## 4 Bolivia Bolivia BO South America 11513100 2019
## 5 Peru Peru PE South America 32510453 2019
## 6 Argentina Argentina AR South America 44938712 2019
## geometry
## 1 MULTIPOLYGON (((117.7036 4....
## 2 MULTIPOLYGON (((117.7036 4....
## 3 MULTIPOLYGON (((-69.51009 -...
## 4 MULTIPOLYGON (((-69.51009 -...
## 5 MULTIPOLYGON (((-69.51009 -...
## 6 MULTIPOLYGON (((-67.1939 -2...## Simple feature collection with 0 features and 6 fields
## Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
## Geodetic CRS: WGS 84
## [1] SOVEREIGNT ADMIN ISO_A2 CONTINENT POP_EST POP_YEAR geometry
## <0 rows> (or 0-length row.names)There are no missing values for the population field (an improvement
with respect to the world dataset).
We will make another version of this dataframe with some of the continent data rearranged: * We will drop the continent of Antarctica, as it adds little value to the histograms. * We will group the countries in the group “Seven seas (open ocean)” into their respective continents.
gmsf_population_by_countries_sel_rearranged <- gmsf_population_by_countries_sel %>%
mutate(CONTINENT = case_when(
ADMIN == "French Southern and Antarctic Lands" ~ "Antarctica", # closest to Antarctica and has 140 citizens
ADMIN == 'Seychelles' ~ 'Africa',
ADMIN == 'Heard Island and McDonald Islands' ~ "Antarctica", # closest to Antarctica and has 0 citizens
ADMIN == 'Saint Helena' ~ 'Africa',
ADMIN == 'Mauritius' ~ 'Africa',
ADMIN == 'British Indian Ocean Territory' ~ 'Asia',
ADMIN == 'Maldives' ~ 'Asia',
ADMIN == 'South Georgia and the Islands' ~ 'South America',
ADMIN == 'Clipperton Island' ~ 'North America', # as Central America is included in North America in this dataset
TRUE ~ CONTINENT
)) %>%
filter(CONTINENT != "Antarctica")
head(gmsf_population_by_countries_sel_rearranged)## Simple feature collection with 6 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -109.4537 ymin: -55.9185 xmax: 140.9776 ymax: 7.35578
## Geodetic CRS: WGS 84
## SOVEREIGNT ADMIN ISO_A2 CONTINENT POP_EST POP_YEAR
## 1 Indonesia Indonesia ID Asia 270625568 2019
## 2 Malaysia Malaysia MY Asia 31949777 2019
## 3 Chile Chile CL South America 18952038 2019
## 4 Bolivia Bolivia BO South America 11513100 2019
## 5 Peru Peru PE South America 32510453 2019
## 6 Argentina Argentina AR South America 44938712 2019
## geometry
## 1 MULTIPOLYGON (((117.7036 4....
## 2 MULTIPOLYGON (((117.7036 4....
## 3 MULTIPOLYGON (((-69.51009 -...
## 4 MULTIPOLYGON (((-69.51009 -...
## 5 MULTIPOLYGON (((-69.51009 -...
## 6 MULTIPOLYGON (((-67.1939 -2...1.4. Ports data
Source of the ports data: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/ports/ Select only the values that may be useful:
gmsf_ports <- st_read("/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/geo_assig2/map_1/ports/ne_10m_ports.shp")## Reading layer `ne_10m_ports' from data source
## `/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/geo_assig2/map_1/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 841.5. Airports data
Source of the airports data: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/airports/
Select only the features that may be useful:
gmsf_airports <- st_read("/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/geo_assig2/map_1/airports/ne_10m_airports.shp")## Reading layer `ne_10m_airports' from data source
## `/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/geo_assig2/map_1/airports/ne_10m_airports.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 893 features and 40 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -175.1356 ymin: -53.78147 xmax: 179.1954 ymax: 78.24672
## Geodetic CRS: WGS 84## [1] "scalerank" "featurecla" "type" "name" "abbrev"
## [6] "location" "gps_code" "iata_code" "wikipedia" "natlscale"
## [11] "comments" "wikidataid" "name_ar" "name_bn" "name_de"
## [16] "name_en" "name_es" "name_fr" "name_el" "name_hi"
## [21] "name_hu" "name_id" "name_it" "name_ja" "name_ko"
## [26] "name_nl" "name_pl" "name_pt" "name_ru" "name_sv"
## [31] "name_tr" "name_vi" "name_zh" "wdid_score" "ne_id"
## [36] "name_fa" "name_he" "name_uk" "name_ur" "name_zht"
## [41] "geometry"2. Maps
2.1. Map of total population by country
gmsf_population_by_countries_sel <- gmsf_population_by_countries_sel %>%
mutate(POP_EST_Millions = POP_EST / 1e6)
ggplot(gmsf_population_by_countries_sel) +
geom_sf(aes(fill=POP_EST_Millions)) + # Fill countries by population estimates
scale_fill_distiller(palette = "Spectral", # Setting color legend
name = "Population (in Millions)") +
labs(title = "World Population by Country",
subtitle = "Population estimates across countries",
caption = "Data Source: Natural Earth (2025)")
We scaled the data to millions for a more clear legend. Note, however, that the population estimates correspond to different years, so the current situation may differ slightly.
2.2. Histogram of country population distribution by continent
For this, we used the rearranged population by country, to avoid having a histogram with just a 6-7 values and one for Antarctica alone.
gmsf_population_by_countries_sel_rearranged <- gmsf_population_by_countries_sel_rearranged %>%
mutate(POP_EST_Millions = POP_EST / 1e6)
ggplot(gmsf_population_by_countries_sel_rearranged, aes(POP_EST_Millions, fill = CONTINENT)) +
geom_histogram(bins = 30, color = "white") +
facet_wrap(~ CONTINENT, scales = "free_x") +
scale_fill_manual(values = c(
"Africa" = "darkgreen",
"Asia" = "darkblue",
"Europe" = "darkred",
"North America" = "darkorange",
"South America" = "darkviolet",
"Oceania" = "darkcyan"
)) +
labs(
title = "Population by Continent (different scales)",
x = "Population Estimate (in Millions)",
y = "Frequency"
) +
theme_minimal() +
theme(
strip.text = element_text(size = 12, face = "bold"),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8),
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10),
legend.position = "none"
)
We also plot it in log scale, to make it more readable.
ggplot(gmsf_population_by_countries_sel_rearranged, aes(POP_EST, fill = CONTINENT)) +
geom_histogram(bins = 30, color = "white") +
facet_wrap(~ CONTINENT, scales = "free_x") +
scale_fill_manual(values = c(
"Africa" = "darkgreen",
"Asia" = "darkblue",
"Europe" = "darkred",
"North America" = "darkorange",
"South America" = "darkviolet",
"Oceania" = "darkcyan"
)) +
scale_x_log10(
breaks = c(100, 1000, 10000, 100000, 1000000, 10000000, 10000000, 1000000000),
labels = c('100', '1k', '10k', '100k', '1M', '10M', '100M', '1B')
) +
labs(
title = "Population Histogram by Continent (Log Scale)",
x = "Log of Population Estimate",
y = "Frequency"
) +
scale_y_continuous(
breaks = seq(0, 13, by = 3)
) +
theme_minimal() +
theme(
strip.text = element_text(size = 12, face = "bold"),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8),
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)
)## Warning in scale_x_log10(breaks = c(100, 1000, 10000, 1e+05, 1e+06, 1e+07, :
## log-10 transformation introduced infinite values.## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).
2.3. Histogram of (country-level) average distances between locations and ports or airports by continent
For creating these histograms, we need to: 1. First, compute the distances of each location to all the airports and ports within the same country (if computationally feasible, if not limit to top 20 locations by country). 2. Second, keep just the shortest distance of each location within a country to an airport and a port. 3. Third, compute the average distances (country-level) of each location to an airport and a port. 4. Fourth, plot the histogram with the average distances by continent.
Preliminary cleaning
First, we filter out those airports which only have a military or spaceport purpose:
## [1] "small" "mid" "mid and military"
## [4] "major and military" "military mid" "military"
## [7] "major" "military major" "spaceport"We remove those airports which just have a military or spaceport purpose, since they won’t be an indicator for how well connected a location is.
gmsf_airports_slice <- gmsf_airports_sel %>%
filter(!(type %in% c("military", "military mid", "military major", "spaceport")))
unique(gmsf_airports_slice$type)## [1] "small" "mid" "mid and military"
## [4] "major and military" "major"Step 1: Compute distances of (populated) locations to infrastructures in the same country
What is the number of calculations that we would have to do if we wanted to find the distances of all locations from all the infrastructures?
print(paste("Number of calculations for airports:", dim(gmsf_airports_slice)[1] * dim(gmsf_population_sel)[1]))## [1] "Number of calculations for airports: 6416908"print(paste("Number of calculations for ports:", dim(gmsf_ports_sel)[1] * dim(gmsf_population_sel)[1]))## [1] "Number of calculations for ports: 7936702"So, in the case of airports, we will have to do 6,416,908 distance calculations. For ports, 7,936,702.
Calculations for ports:
The resulting matrix has dim(gmsf_population_sel)[1]
rows and dim(gmsf_ports_sel)[1] columns.
Each element [i, j] in the matrix represents the distance between the i-th geometry in the first object and the j-th geometry in the second object.
Calculations for airports:
Step 2: Keep the shortest distance of each location with each type of infrastructure
Since the resulting matrices have
dim(gmsf_population_sel)[1] rows and
dim(gmsf_ports_sel)[1] columns, we can select the minimum
distance of each location with each infrastructure by selecting the
minimum value of each row.
## [1] 7342 874## [1] 7342 1081Now, we compute the minimum distance for each row (location) with each type of infrastructure, and save it as a vector:
min_airports_distances <- apply(
X = airports_distance, MARGIN = 1, FUN = min, na.rm = TRUE
)
min_ports_distances <- apply(
X = ports_distance, MARGIN = 1, FUN = min, na.rm = TRUE
)We convert the vectors to km for making the units more readable:
min_airports_distances_km <- min_airports_distances / 1000
min_ports_distances_km <- min_ports_distances / 1000
head(min_airports_distances_km)## [1] 55.39118 42.56526 89.21916 41.09979 84.29933 347.79542Below, we bind the minimum distance vectors of each type of
infrastructure as new columns of the gmsf_population_sel
sf.
Since the order of the rows in gmsf_population_sel
remains unchanged during the distance matrix calculation, we can
directly bind the minimum distance vector as a new column.
gmsf_population_sel$min_distance_airport_km <- min_airports_distances_km
gmsf_population_sel$min_distance_port_km <- min_ports_distances_km
head(gmsf_population_sel)## Simple feature collection with 6 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58.304 ymin: -34.538 xmax: 0.7890036 ymax: 9.261
## Geodetic CRS: WGS 84
## FEATURECLA NAME NAMEASCII SOV0NAME
## 1 Admin-1 capital Colonia del Sacramento Colonia del Sacramento Uruguay
## 2 Admin-1 capital Trinidad Trinidad Uruguay
## 3 Admin-1 capital Fray Bentos Fray Bentos Uruguay
## 4 Admin-1 capital Canelones Canelones Uruguay
## 5 Admin-1 capital Florida Florida Uruguay
## 6 Admin-1 capital Bassar Bassar Togo
## ADM0NAME ISO_A2 POP_MAX POP_MIN geometry
## 1 Uruguay UY 21714 21714 POINT (-57.83612 -34.46979)
## 2 Uruguay UY 21093 21093 POINT (-56.901 -33.544)
## 3 Uruguay UY 23279 23279 POINT (-58.304 -33.139)
## 4 Uruguay UY 19698 19698 POINT (-56.284 -34.538)
## 5 Uruguay UY 32234 32234 POINT (-56.215 -34.099)
## 6 Togo TG 61845 61845 POINT (0.7890036 9.261)
## min_distance_airport_km min_distance_port_km
## 1 55.39118 0.5927656
## 2 42.56526 134.7817764
## 3 89.21916 124.9193273
## 4 41.09979 40.9677677
## 5 84.29933 89.1342812
## 6 347.79542 351.4064671Step 3: Compute average distances (country-level) of each location to each type of infrastructure
Now, the idea is to group by the sovereign country name and compute the average distance of each location with each type of infrastructure, ports and airports:
# Group by country (using the ISO_A2 country code)
avg_country_dist_airport <- gmsf_population_sel %>%
group_by(ISO_A2) %>%
summarise(avg_dist_airport_km = mean(min_distance_airport_km, na.rm = T))
avg_country_dist_port <- gmsf_population_sel %>%
group_by(ISO_A2) %>%
summarise(avg_dist_port_km = mean(min_distance_port_km, na.rm = T))## Simple feature collection with 6 features and 2 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -61.85003 ymin: 8.433297 xmax: 71.15 ymax: 42.66671
## Geodetic CRS: WGS 84
## # A tibble: 6 × 3
## ISO_A2 avg_dist_airport_km geometry
## <chr> <dbl> <GEOMETRY [°]>
## 1 -99 212. MULTIPOINT ((20.31074 42.66033), (20.75009 42.2293…
## 2 AD 125. POINT (1.526594 42.51075)
## 3 AE 48.5 MULTIPOINT ((54.36659 24.46668), (55.01074 24.9762…
## 4 AF 219. MULTIPOINT ((70.57925 37.12976), (68.87253 36.7279…
## 5 AG 6.61 POINT (-61.85003 17.11804)
## 6 AL 63.5 MULTIPOINT ((19.49823 40.47736), (19.51885 42.0684…## Simple feature collection with 6 features and 2 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -61.85003 ymin: 8.433297 xmax: 71.15 ymax: 42.66671
## Geodetic CRS: WGS 84
## # A tibble: 6 × 3
## ISO_A2 avg_dist_port_km geometry
## <chr> <dbl> <GEOMETRY [°]>
## 1 -99 154. MULTIPOINT ((20.31074 42.66033), (20.75009 42.22932),…
## 2 AD 139. POINT (1.526594 42.51075)
## 3 AE 38.1 MULTIPOINT ((54.36659 24.46668), (55.01074 24.97624),…
## 4 AF 1097. MULTIPOINT ((70.57925 37.12976), (68.87253 36.72795),…
## 5 AG 0.579 POINT (-61.85003 17.11804)
## 6 AL 74.5 MULTIPOINT ((19.49823 40.47736), (19.51885 42.06845),…Now, we drop the geometries from the data frames containing the
average distances and merge the results with the
gmsf_population_by_countries_sel data frame, which has the
same names for the sovereign countries and also contains the country
polygons:
# First, we drop the geometries
avg_country_dist_airport <- st_drop_geometry(avg_country_dist_airport)
avg_country_dist_port <- st_drop_geometry(avg_country_dist_port)
# Second, we do a left join on the distance data frames, in order to
# have the polygons of the world data frame
gmsf_avg_distances <- left_join(
x = gmsf_population_by_countries_sel_rearranged,
y = avg_country_dist_airport,
by = "ISO_A2")
gmsf_avg_distances <- left_join(
x = gmsf_avg_distances,
y = avg_country_dist_port,
by = "ISO_A2")
head(gmsf_avg_distances)## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -109.4537 ymin: -55.9185 xmax: 140.9776 ymax: 7.35578
## Geodetic CRS: WGS 84
## SOVEREIGNT ADMIN ISO_A2 CONTINENT POP_EST POP_YEAR POP_EST_Millions
## 1 Indonesia Indonesia ID Asia 270625568 2019 270.62557
## 2 Malaysia Malaysia MY Asia 31949777 2019 31.94978
## 3 Chile Chile CL South America 18952038 2019 18.95204
## 4 Bolivia Bolivia BO South America 11513100 2019 11.51310
## 5 Peru Peru PE South America 32510453 2019 32.51045
## 6 Argentina Argentina AR South America 44938712 2019 44.93871
## avg_dist_airport_km avg_dist_port_km geometry
## 1 251.8960 178.43623 MULTIPOLYGON (((117.7036 4....
## 2 122.0240 72.57155 MULTIPOLYGON (((117.7036 4....
## 3 161.4024 105.07845 MULTIPOLYGON (((-69.51009 -...
## 4 200.5321 613.25985 MULTIPOLYGON (((-69.51009 -...
## 5 287.5093 242.84726 MULTIPOLYGON (((-69.51009 -...
## 6 224.4619 349.97085 MULTIPOLYGON (((-67.1939 -2...Step 4: Plotting histograms by type of infrastructure
Note that not all locations in the world are included in the dataset that has been used, so the histograms may not accurately represent reality.
ggplot(gmsf_avg_distances, aes(x = avg_dist_airport_km, fill = CONTINENT)) +
geom_histogram(bins = 30, color = "white", na.rm = TRUE) +
facet_wrap(~ CONTINENT, scales = "free_x") +
scale_fill_manual(values = c(
"Africa" = "darkgreen",
"Asia" = "darkblue",
"Europe" = "darkred",
"North America" = "darkorange",
"South America" = "darkviolet",
"Oceania" = "darkcyan"
)) +
labs(
title = "Country-level Average Distances to Closest AIRPORT",
x = "Average Distance (km)",
y = "Frequency"
) +
theme_minimal() +
theme(
strip.text = element_text(size = 12, face = "bold"),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8),
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10),
legend.position = "none"
)
ggplot(gmsf_avg_distances, aes(avg_dist_port_km, fill = CONTINENT)) +
geom_histogram(bins = 30, color = "white", na.rm = TRUE) +
facet_wrap(~ CONTINENT, scales = "fixed") +
labs(title = "Country-level Average Distances to Closest PORT",
x = "Average Distance (km)",
y = "Frequency") +
scale_fill_manual(values = c(
"Africa" = "darkgreen",
"Asia" = "darkblue",
"Europe" = "darkred",
"North America" = "darkorange",
"South America" = "darkviolet",
"Oceania" = "darkcyan"
)) +
theme_minimal() +
theme(
strip.text = element_text(size = 12, face = "bold"),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8),
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10),
legend.position = "none"
)
Additional: maps of country-level average distances to the closest ports and airports
ggplot(gmsf_avg_distances) +
geom_sf(aes(fill=avg_dist_airport_km)) + # Fill countries by avg distances
scale_fill_distiller(palette = "Spectral", # Setting color legend
name = "Average Distance to Airport") +
labs(title = "Average Distance to Closest Airport",
subtitle = "Country-level average distances to the closest airport, in km (using the arithmetic mean for populated locations)",
caption = "Data Source: Natural Earth (2025)",
)
ggplot(gmsf_avg_distances) +
geom_sf(aes(fill=avg_dist_port_km)) +
scale_fill_distiller(palette = "Spectral",
name = "Average Distance to Port") +
labs(title = "Average Distance to Closest Port",
subtitle = "Country-level average distances to the closest port, in km\n(using the arithmetic mean for populated locations)",
caption = "Data Source: Natural Earth (2025)")
\(\color{darkblue}{\text{Part 2}}\)
1. Introduction
This assignment focuses on analyzing market locations in sub-Saharan Africa and their relationship with infrastructure, including roads, airports, and coastlines. We will use geospatial techniques to compute distances and visualize key trends in market prices.
2. Load the location of markets
We start by loading the market location data, which includes latitude and longitude coordinates. The dataset was found in the replication package published by the author.
cwd <- '/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2'
file_name <- 'data/MktCoords.xlsx'
file_path <- file.path(cwd, file_name)
print(file_path)## [1] "/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/data/MktCoords.xlsx"## [1] "/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2"3. Country Boundaries
We have applied several filters to the world dataset: - We noticed that the paper only looks at Sub-Saharan Africa, so we excluded Northern African countries. - We excluded Madagascar, because it’s not present in the paper’s research. - The Republic of Sudan is attributed to Northern Africa in the world dataset, however, in the paper it’s considered a part of sub-Saharan Africa. So we decided to include it to be able to replicate the paper’s findings more closely.
4. Roads
Initially we downloaded the roads network from Natural Earth, however, that dataset did not contain information about countries (only continents). So it was impossible to apply those filters that we identified in the previous step (exclude Northern Africa and Madagascar, but include Sudan). So we found an alternative dataset from ArcGis, which did contain country labels.
We only included Primary roads, because otherwise the entire map was covered with a fine mesh of roads, which were adding too much noise for our later calculations of distances.
#Load the roads data from ArcGis (https://www.arcgis.com/home/item.html?id=ba1cf90a739f41f4b91b26441929918a&view=list&sortOrder=desc&sortField=defaultFSOrder#overview)
roads <- st_read('/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/data/AFR_Infra_Transport_Road.shp/AFR_Infra_Transport_Road.shp')## Reading layer `AFR_Infra_Transport_Road' from data source
## `/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/data/AFR_Infra_Transport_Road.shp/AFR_Infra_Transport_Road.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 360852 features and 11 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -17.51197 ymin: -34.82845 xmax: 50.28193 ymax: 37.30094
## Geodetic CRS: WGS 84#Filter road types to keep only primary roads, otherwise it looks too messy
roads_type <- c("Road (Primary)")
main_roads <- roads %>%
filter(FeatureTyp %in% roads_type)
# Exclude Northern African countries
ssa_exceptions <- c("Algeria", "Egypt", "Morocco and Western Sahara", "Madagascar", "Tunisia", "Libya")
roads_ssa <- main_roads %>%
filter(!Country %in% ssa_exceptions)
# Even when filtered, there are some roads that transcend the boundaries of our region of interest, so we use st_within function.
roads_within <- st_within(roads_ssa, africa)
# Convert the sparse matrix to a list of vectors
roads_within_list <- lapply(roads_within, function(x) if (length(x) == 0) NA_integer_ else x)
# Filter the linestrings
roads_within_ssa <- roads_ssa[!is.na(roads_within_list), ]5. Airports
Again, we only kept big airports and excluded small airports and helipads, because usually they are not used for cargo transportation.
#Load airports data (https://ourairports.com/markets/)
airports_csv <- read_csv('/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/data/airports.csv')
#Transform to geometry
airports <- st_as_sf(airports_csv, coords = c("longitude_deg", "latitude_deg"), crs = 4326)
# Filter airports that fall within the polygons of africa dataset
points_within <- st_within(airports, africa)
# Convert the sparse matrix to a list of vectors
points_within_list <- lapply(points_within, function(x) if (length(x) == 0) NA_integer_ else x)
#Select points that are within the polygon
airports_africa <- airports[which(!is.na(points_within_list)), ]
# Filter out small, medium and closed airports
big_airports <- c("large_airport")
big_airports_africa <- airports_africa %>%
filter(type %in% big_airports)6. Prices
We then loaded the price data and computed the average price across all available columns, for each crop for each market.
# Price markets
prices <- read_excel("/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/data/PriceMaster4GAMS.xlsx")
# Compute the average price across columns 1 to 46
prices <- prices %>%
rowwise() %>%
mutate(price_avg = mean(c_across(`1`:`46`), na.rm = TRUE)) %>%
select(mktcode, country, market, crop, price_avg) # Keep relevant columns
# Merge based on market code
# Ensure column names are consistent before merging
colnames(prices)## [1] "mktcode" "country" "market" "crop" "price_avg"## [1] "ctrycode" "mktcode" "market" "geometry"7. Coastline
The coastline data came from Natural Earth.
# Coastline data
coastline <- st_read("/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/data/ne_10m_coastline/ne_10m_coastline.shp")## Reading layer `ne_10m_coastline' from data source
## `/Users/newmac/Documents/DSDM/Term_2/03_GIS/Assig_2/data/ne_10m_coastline/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 848. Distances from markets to objects of infrastructure
In the following block of code we use the “st_distance” function to calculate the distances to various objects of infrastructure: coastline, roads and airports.
# Compute distance to the nearest coastline
dist_matrix_coast <- st_distance(markets, africa_coastline, by_element = FALSE)
markets <- markets %>%
mutate(dist_coast = apply(dist_matrix_coast, 1, min, na.rm = TRUE))
# Compute distance to the nearest road
dist_matrix_roads <- st_distance(markets, roads_within_ssa, by_element = FALSE)
markets <- markets %>%
mutate(dist_road = apply(dist_matrix_roads, 1, min, na.rm = TRUE))
# Compute distance to the nearest airport
dist_matrix_airports <- st_distance(markets, big_airports_africa, by_element = FALSE)
markets <- markets %>%
mutate(dist_airport = apply(dist_matrix_airports, 1, min, na.rm = TRUE)) 9. Visualizations
Simple ggplot with 5 layers to plot all our data on the same map and replicate the figure we found in the paper.
#Plot the markets
ggplot() +
# Layer 1: Countries with boundaries
geom_sf(data = africa, fill = "khaki") +
# Layer 2: Coastline
geom_sf(data = africa_coastline, aes(color = "Coastline"), size = 0.6, alpha = 0.7) +
# Layer 3: Roads
geom_sf(data = roads_within_ssa, size = 0.1, aes(color = "Roads"), alpha = 0.8) +
# Layer 4: Market locations
geom_sf(data = markets, aes(color = "Markets"), size = 0.8) +
# Layer 5: Airports
geom_sf(data = big_airports_africa, aes(color = "Airports"), size = 2, shape = 17) +
labs(title = "Market Locations and Infrastructure Across Africa",
x = "Longitude", y = "Latitude",
color = "Infrastructure") +
scale_color_manual(values = c("Markets" = "black",
"Coastline" = "blue",
"Roads" = "orange",
"Airports" = "purple")) +
theme_minimal() +
theme(
axis.title = element_blank(),
axis.text = element_blank(),
legend.title = element_text(face = "bold"),
plot.title = element_text(hjust = 0.5, face = "bold", size = 14)
)
10. Scatter plots of log(distance) vs crop prices
We decided to build two sets of scatter plots for each type of infrastructure (coastline, roads and airports): separate faceted scatter plots for each crop (using facet_wrap function), and one general scatter plot combining all the crops (but still color-coded).
markets <- markets %>%
filter(!is.na(price_avg) & price_avg > 0) %>%
mutate(
log_dist_coast = log1p(dist_coast),
log_dist_road = log1p(dist_road),
log_dist_airport = log1p(dist_airport),
log_price = log1p(price_avg)
)# Scatter plot of Price vs Distance to Coast by Crop
ggplot(markets, aes(x = log_dist_coast, y = price_avg, color = crop)) +
geom_point(alpha = 0.6) +
facet_wrap(~crop, scales = "free") +
labs(title = "Average Price vs Log Distance to Coast by Crop",
x = "Log Distance to Coast (meters)",
y = "Average Price") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
# Scatter plot of Price vs Distance to Coast All Crops Combined
ggplot(markets, aes(x = log_dist_coast, y = price_avg, color = crop)) +
geom_point(alpha = 0.6) +
labs(title = "Average Price vs Log Distance to Coast by Crop",
x = "Log Distance to Coast (meters)",
y = "Average Price") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
# Scatter plot of Price vs Distance to Nearest Road by Crop
ggplot(markets, aes(x = log_dist_road, y = price_avg, color = crop)) +
geom_point(alpha = 0.6) +
facet_wrap(~crop, scales = "free") +
labs(title = "Average Price vs Log Distance to Road by Crop",
x = "Log Distance to Nearest Road (meters)",
y = "Average Price") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
# Scatter plot of Price vs Distance to Nearest Road All Crops Combined
ggplot(markets, aes(x = log_dist_road, y = price_avg, color = crop)) +
geom_point(alpha = 0.6) +
labs(title = "Average Price vs Log Distance to Road by Crop",
x = "Log Distance to Nearest Road (meters)",
y = "Average Price") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
# Scatter plot of Price vs Distance to Nearest Airport by Crop
ggplot(markets, aes(x = log_dist_airport, y = price_avg, color = crop)) +
geom_point(alpha = 0.6) +
facet_wrap(~crop, scales = "free") +
labs(title = "Average Price vs Log Distance to Airport by Crop",
x = "Log Distance to Nearest Airport (meters)",
y = "Average Price") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
# Scatter plot of Price vs Distance to Nearest Airport All Crops Combined
ggplot(markets, aes(x = log_dist_airport, y = price_avg, color = crop)) +
geom_point(alpha = 0.6) +
#facet_wrap(~crop, scales = "free") +
labs(title = "Average Price vs Log Distance to Airport by Crop",
x = "Log Distance to Nearest Airport (meters)",
y = "Average Price") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
11. Conclusions
We find that crop prices generally rise with distance from the coast, supporting the idea that trade costs increase prices form the paper by Porteous (2019). However, some coastal markets still show high prices, this could be suggesting that there other factors at play for example the quality of the infrastructure. Airports, as noted in the paper, have seem to have little impact.
These finding are consistent with the author’s conclusions, who writes: “I found that lower agricultural trade costs would have led to a large drop in grain prices, agricultural revenues, and expenditure on grains in sub-Saharan Africa during the study period, with an overall welfare gain equivalent to 2.17 percent of GDP. There was significant variation in these effects, with some markets experiencing increases in prices, revenue, and welfare, and others experiencing welfare losses due to terms-of-trade effects.”
Another explanation for this phenomenon could be that markets closer to ports and airports have a higher ratio of imported goods, which are commonly more expensive due to tariffs and transportation costs. It can also be that markets closer to transportation hubs experience higher demand, plus they are concentrated in urban areas where prices are higher in general.
References
Porteous, O., 2019. High trade costs and their consequences: An estimated dynamic model of African agricultural storage and trade. American Economic Journal: Applied Economics, 11(4), pp.327-66.
Data Sources
Market and price data: Paper replication package.
Road data: ArcGis. https://www.arcgis.com/home/item.html?id=ba1cf90a739f41f4b91b26441929918a&view=list&sortOrder=desc&sortField=defaultFSOrder#overview)
Airports data: OurAirports. https://ourairports.com/data/
Coastline data: Natural Earth. https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-coastline/