Making maps in R using ggplot and sf
Setup
# Activating packages
suppressWarnings(suppressMessages({
library(tidyverse) # Modern data science workflow
library(stringr) # Functions designed to work with strings
library(sf) # Encodes spatial vector data
library(RColorBrewer) # ColorBrewer Palettes
library(leaflet) # Interactive web maps
library(ggplot2) # Modern visualization tools
library(ggthemes) # For extra themes, scales and geoms for 'ggplot2'
library(patchwork) # make it simple to combine separate ggplots
library(plotly) # Create interactive web graphics from 'ggplot2'
library(ggrepel) ## For displaying labels on ggplot2 object
}))
# Adjustments
knitr::opts_chunk$set(fig.align = 'center')
options(warn = -1)
options(scipen = 10000)
options(repr.plot.width = 6, repr.plot.height = 4)
Import data
df <- read.csv("https://raw.githubusercontent.com/haginta/regional-convergence-indonesia-satellite-gdp/main/java118_.csv")
glimpse(df)
mymap <- read_sf("https://github.com/haginta/regional-convergence-indonesia-satellite-gdp/raw/main/java_118.gpkg")
glimpse(mymap)
st_crs(mymap)
Transform data
df[,4:7] <- log(df[4:7])
head(df)
df <- df %>%
mutate(sat_g = sat2020 - sat2019,
gdp_g = gdp2020 - gdp2019)
head(df)
Merge datasets
df_map <- merge(mymap, df, by=c("districtID", "district", "province"))
df_map <- st_make_valid(df_map)
glimpse(df_map)
Maps
Simple boundaries
df_map %>%
ggplot() +
geom_sf()
Labeled boundaries
#Create dataframe of selected cities
city <- c("Kota Jakarta Pusat", "Kota Yogyakarta", "Kota Surabaya")
h_city <- df_map %>% filter(district %in% city)
#Add coordinates
coordinates <- as.data.frame(st_coordinates(st_centroid(h_city)))
h_city$X <- coordinates$X
h_city$Y <- coordinates$Y
glimpse(h_city)
ggplot() +
geom_sf(data = h_city, aes(fill = district), lwd = 0) +
geom_sf(data = df_map, aes(fill = NA, alpha=0), show.legend = F) +
geom_sf(data = h_city, aes(fill = district, alpha=0), show.legend = F, color="red") +
geom_label_repel(data = h_city, aes(x = X, y = Y, label = district), size=2, color="blue", point.padding = NA) +
theme(legend.position="None") +
labs(x = "Longitude", y="Latitude")
Choropleth map
ggplot(df_map) +
geom_sf(aes(fill=sat2019), color = "white", lwd = 0.05) +
scale_fill_viridis_c(option = "magma", name = "Radiance (log)") +
theme_map() +
theme(legend.direction="horizontal") +
labs(title = "Districts of Java seen from space at night",
subtitle = "Nighttime light (NTL) intensity across 118 districts in Java, Indonesia, 2019",
caption = "Geometries: GADM \n Data: Earth Observation Group")
ggplot() +
geom_sf(data = h_city, aes(fill = sat2019), color = "white", lwd = 0.05) +
geom_sf(data = df_map, aes(fill = sat2019), color = "white", lwd = 0.05, show.legend = T) +
geom_sf(data = h_city, aes(fill = sat2019), show.legend = F, color="red") +
scale_fill_viridis_c(option = "magma", name = "Radiance (log)") +
theme_map() +
geom_label_repel(data = h_city, aes(x = X, y = Y, label = district), size=2, color="blue", point.padding = NA) +
theme(legend.direction="horizontal", legend.position=c(0.05, -0.12))+
labs(x = "Longitude", y="Latitude")
Comparative maps
p1 <- ggplot(df_map) +
geom_sf(aes(fill=sat2020), color = "white", lwd = 0.05) +
scale_fill_viridis_c(option = "magma", name = "Radiance (log)") +
theme_map() +
theme(legend.direction="horizontal", legend.position=c(-0.0, -0.15))+
labs(subtitle = "(a) Spatial distribution of NTL per capita")
p2 <- ggplot(data = df_map) +
geom_sf(aes(fill=gdp2020), color = "white", lwd = 0.05) +
scale_fill_viridis_c(option = "magma", name = "GDP (log)") +
theme_map() +
theme(legend.direction="horizontal", legend.position=c(-0.0, -0.15))+
labs(subtitle = "(b) Spatial distribution of GDP per capita")
p1 / p2
Interactive
p <- ggplot(df_map) + geom_sf(aes(fill=gdp2019, label=district), lwd = 0.05) +
scale_fill_viridis_c(option = "magma", name = "GDP (log)") + theme_bw() +
labs(x = "Longitude", y="Latitude", title = "GDP across 118 districts in Java, Indonesia, 2019 (constant price 2010)")
p_inter <- ggplotly(p)
p_inter