-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path7.bite.sized.vis.1.ind_bl_pol_dist.Rmd
203 lines (160 loc) · 6.46 KB
/
7.bite.sized.vis.1.ind_bl_pol_dist.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
---
title: "7.bite.sized.vis.1.ind_bl_pol_dist"
author: "Aarsh Batra"
date: "`r Sys.Date()`"
output: html_document
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
## load libraries and background helper functions
```{r load_lib}
library(tidyverse)
library(jsonlite)
library(here)
library(sf)
library(ggplot2)
library(geojson)
library(devtools)
library(geojsonio)
library(geojsonsf)
library(roxygen2)
library(magrittr)
library(stringr)
library(data.table)
library(viridis)
library(ggthemes)
library(hrbrthemes)
library(RColorBrewer)
library(haven)
library(janitor)
library(lwgeom)
library(rvest)
library(dplyr)
library(ggrepel)
library(ggsflabel)
library(fuzzyjoin)
library(tidytext)
library(tm)
library(wordcloud2)
library(recipes)
library(forcats)
library(forcats)
library(geofacet)
library(pak)
library(camcorder)
library(ggbeeswarm)
library(ggtext)
library(ggimage)
library(tidyverse)
library(ggbeeswarm)
library(cowplot)
library(ggtext)
library(scales)
library(grid)
library(png)
library(here)
library(showtext)
library(sysfonts)
library(ggimage)
library(Cairo)
library(extrafont)
library(MetBrewer)
library(rnaturalearth)
library(scico)
# Negation of the "%in%" operator function
`%notin%` <- Negate(`%in%`)
# load SHRUG block level shp file
ind_shr_subdist_shp <- st_read(paste0(here(), "/4.ind.bl.lev.pol.dist/raw.data/shrug_master_subdistrict_shp.gpkg"))
ind_shr_subdist_shp <- ind_shr_subdist_shp %>%
mutate_if(is.character, tolower)
# load india shp
ind_country_shp <- st_read(paste0(here(), "/4.ind.bl.lev.pol.dist/raw.data/india_country_shrug.gpkg"))
```
## India choropleth 2022 satellite derived PM2.5 pollution
```{r}
#### load raw data
ind_bl_lev_pol_wide <- read_csv(paste0(here(), "/4.ind.bl.lev.pol.dist/raw.data/ind_block_level_pop_weighted_pol_1998_2022.csv"))
#### reshape and keep only year 2022
ind_bl_lev_pol_long_2022 <- ind_bl_lev_pol_wide %>%
pivot_longer(cols = starts_with("avg_pm"), names_to = "pol_year", values_to = "ann_avg_pm2.5") %>%
mutate(pol_year = as.numeric(str_remove(str_extract(pol_year, "_\\d+"), "_")))
#### create a 2022 only wide version
ind_bl_lev_pol_wide_2022 <- ind_bl_lev_pol_wide %>%
select(pc11_state_id:subdistrict_population, avg_pm2.5_2022)
#### join 2022 wide with shrug block shp
ind_bl_lev_pol_wide_2022_shp <- ind_bl_lev_pol_wide_2022 %>%
left_join(ind_shr_subdist_shp, by = c("pc11_state_id", "pc11_district_id",
"pc11_subdistrict_id", "state_name",
"district_name", "subdistrict_name")) %>%
select(-geom, geom) %>%
st_as_sf() %>%
mutate(pollution_level = cut(
avg_pm2.5_2022,
breaks = c(0, 5, 20, 40, 80, 120), # Set the breaks for pollution levels
labels = c("Complying with WHO", "5 to 20", "20 to 40", "40 to 80", "80+"), # Adjust labels to match the number of intervals
include.lowest = TRUE,
right = FALSE # Open intervals on the right side, so a value of 5 falls in the second category
)) %>%
select(-geom, geom)
#### collapsing by pollution level
ind_bl_lev_pol_wide_2022_shp_summary <- ind_bl_lev_pol_wide_2022_shp %>%
group_by(pollution_level) %>%
summarise(geom = st_union(geom)) %>%
st_simplify()
#### Define color palette (use colors that get more concerning as pollution increases)
colors_grey <- c("#FFFFFF", "#707070", "#3C3C3C", "#2B2B2B", "black")
#### load fonts
font_import()
loadfonts()
windowsFonts()
f1 <- "Baskerville Old Face"
font_add_google(f1)
showtext_auto(FALSE)
#### Plot the map with a stepped color scale
plt <- ggplot(ind_bl_lev_pol_wide_2022_shp) +
# Add country boundary with transparent fill and thin black outline
geom_sf(data = ind_country_shp, fill = "transparent", color = "#CCCCCC", linewidth = 0.2) +
geom_sf(data = ind_bl_lev_pol_wide_2022_shp_summary, fill = "transparent", color = "#CCCCCC", linewidth = 0.2) +
# Map PM2.5 data with transparent borders to avoid clutter
geom_sf(mapping = aes(fill = avg_pm2.5_2022), color = "transparent") +
# Define fill scale with grayscale color palette and custom breaks/labels
scale_fill_stepsn(
name = expression(PM[2.5] ~ "(µg/m³)"),
colors = colors_grey,
breaks = c(0, 5, 20, 40, 80), # Define thresholds
labels = c("0", "5 (WHO Guideline)", "20", "40 (National Standard)", "80"),
limits = c(0, 110),
oob = scales::squish
) +
# Ensure map coordinates are fitted properly
coord_sf(expand = FALSE) +
# Add titles, subtitles, and captions
labs(
title = "When the Air Chokes: Can India Break Free from Pollution’s Grip?",
subtitle = str_wrap(
"2022 satellite-derived PM2.5 data shows all 5,969 Indian blocks, home to over 1.3 billion people, exceed the WHO annual avg PM2.5 guideline. Nearly half of these blocks, housing over 900 million people (65% of India’s population), exceed the National Standard, with 20% of this population (approx. 180 million) exposed to levels over twice that limit, more than 16 times the WHO guideline",
width = 90
),
caption = expression("Graphic · Aarsh Batra | Pol Data · ACAG, U Wash. 2022" ~ PM[2.5] ~ "Sat. Data | Shp file · SHRUG, DDL | Processed · github.com/AarshBatra/biteSizedAQ")
) +
# Theme adjustments for map and text elements
theme_void() + # Use a blank theme to highlight map features only
theme(
plot.title = element_text(size = 34, family = f1, hjust = 0, color = "black", face = "bold", margin = margin(b = 0.3, t = 2.2, l = 1, r = 0.8, unit = "cm")), # Dark gray title
plot.subtitle = element_text(size = 22, family = f1, hjust = 0, color = "#202020", margin = margin(b = 1, l = 1, r = 1, unit = "cm"), face = "bold.italic"), # Medium gray subtitle
plot.caption = element_text(family = f1, size = 15, color = "#ffffff", margin = margin(t = 2, b = 2, l = 1, r = 1, unit = "cm"), hjust = 0), # Lighter gray caption
legend.text = element_text(family = f1, size = 20, color = "#CCCCCC"), # Mid-gray legend text
legend.title = element_text(family = f1, size = 20, color = "#CCCCCC", hjust = 0), # Mid-gray legend text
legend.position.inside = c(1.8, 0.7),
legend.position = c(1.4, 0.7),
legend.key.height = unit(1.5, "cm"), # Adjust height for better visibility
legend.spacing.y = unit(0.5, "cm")
# plot.margin = margin(t = 10, r = 20, b = 10, l = 10)
)
plt1 <- ggbackground(plt, "./4.ind.bl.lev.pol.dist/images/light_smoke_2.png")
ggsave("./plt1.pdf", plt1, width = 16, height = 10, dpi = 320)
```