Assignment 7: Air Polution
Marissa Beukelman
4/10/2020
# Loading libraries
library(dplyr)
library(rgdal)
library(countrycode)
library(ggmap)
library(readxl)
library(leaflet)
library(DT)
library(taucharts)
library(sp)1. An interactive table of the top 100 cities in the world with the city with the highest PM 10 levels on top.
# Reading the file PM10 Data from Sheet 4
cities = read_excel("OAP_database.xls", sheet = 4)# Elimination rows and retaining only 3 columns
cities <- cities[-c(1:2),c(3,2,4)]# Renaming the 3 columns
names(cities) <- c("City","Country","PMLevel")# Viewing the structure of the cities
str(cities)## tibble [1,099 × 3] (S3: tbl_df/tbl/data.frame)
## $ City : chr [1:1099] "Alger" "Buenos Aires" "Adelaide" "Brisbane" ...
## $ Country: chr [1:1099] "Algeria" "Argentina" "Australia" "Australia" ...
## $ PMLevel: chr [1:1099] "42" "38" "13.5" "18.166666666666668" ...The PMLevel should be considered a number, not a character.
# Changing PMLevel to numeric
cities$PMLevel <- round(as.numeric(cities$PMLevel),2)# Previewing the data of cities
head(cities)## # A tibble: 6 x 3
## City Country PMLevel
## <chr> <chr> <dbl>
## 1 Alger Algeria 42
## 2 Buenos Aires Argentina 38
## 3 Adelaide Australia 13.5
## 4 Brisbane Australia 18.2
## 5 Bunbury Australia 13.8
## 6 Canberra Australia 10.3# Worst 100 cities with high PM levels in an interactive graph
cities %>%
arrange(-PMLevel) %>%
top_n(100) %>%
datatable(cities, rownames = FALSE)2. An interactive bar chart that shows the countries and the number of top 100 cities that are present in each country. The country with the highest number of cities should appear first.
# Counting the number of cities in each country
tmp <- cities %>%
arrange(-PMLevel) %>%
top_n(100) %>%
group_by(Country) %>%
summarise(number_of_cities = length(Country)) %>%
arrange(-number_of_cities)# Creating an interactive graph of top 100 cities by country
tmp$Country <- forcats:: fct_inorder(tmp$Country)
tauchart(tmp) %>%
tau_bar("number_of_cities","Country", horizontal = "TRUE") %>%
tau_legend() %>%
tau_tooltip()3. Plot the top 100 cities in a static world map. (Size the points as a function of the amount of PM 10 present in each city.)
NOTE: I did the Google API, but I will not be paying for the service in order for this assignment to work. I do not have the money to do so at the current time. I did all the coding and connected to a key, but it will not show the dots.
cities100 <- cities %>%
arrange(-PMLevel) %>%
top_n(100)# Creating a new variable CityCountry that combines City and Country separated by a comma
cities100$CityCountry <- paste(cities100$City, cities100$Country, sep=",")# Previewing the cities100
head(cities100)## # A tibble: 6 x 4
## City Country PMLevel CityCountry
## <chr> <chr> <dbl> <chr>
## 1 Ahwaz Iran (Islamic Republic … 372 Ahwaz,Iran (Islamic Republic of)
## 2 Ulaanbaatar Mongolia 279 Ulaanbaatar,Mongolia
## 3 Sanandaj Iran (Islamic Republic … 254 Sanandaj,Iran (Islamic Republic …
## 4 Ludhiana India 251 Ludhiana,India
## 5 Quetta Pakistan 251. Quetta,Pakistan
## 6 Kermanshah Iran (Islamic Republic … 229 Kermanshah,Iran (Islamic Republi…locs <- geocode(cities100$CityCountry)cities100$lat <-locs$lat
cities100$lon <-locs$lonworldmap <- borders("world", fill="lightgray", colour = "white")
worldmap <- ggplot() + worldmap# Adding the cities
worldmap <- worldmap + geom_point(data = cities100, aes(x=lon, y=lat, size=PMLevel),alpha=.4,color="red")
worldmap + theme_void()
4. Plot the top 100 cities in an interactive world map that shows a popup when a city is clicked. The popup should show information on the name of the city, the country of the city, and the annual mean PM 10 level in each city.
cities100$popup <- paste("<table><tr><td>City:", cities100$City, "<br>Country:", cities100$Country, "<br>Annual Mean PM10 Level:", cities100$PMLevel, "</td></tr></table>")leaflet(cities100) %>%
addProviderTiles("CartoDB.Positron") %>%
setView(0, 0, zoom = 2) %>%
addCircles(stroke = FALSE, fillOpacity = .5, color = "red", radius = ~PMLevel*1000, popup = ~popup)5. Using the data for different countries from sheet 5, create an interactive choropleth of the world which colors different countries based on the PM 10 levels. Yet again, the popup when a country is clicked should provide information on the name of the country and the PM10 levels.
# # Reading the file PM10 Data from Sheet 5
countries <- read_excel("OAP_database.xls", sheet = 5)# Elimination rows and retaining 3 columns
countries <- countries[-c(1:2), c(2:3)]# Renaming columns
names(countries) <- c("Country", "PMLevel")# Viewing the structure of countries
str(countries)## tibble [91 × 2] (S3: tbl_df/tbl/data.frame)
## $ Country: chr [1:91] "Estonia" "Mauritius" "Australia" "New Zealand" ...
## $ PMLevel: chr [1:91] "11.132999999999999" "11.65" "13.179790073011779" "15" ...# Changing PMLevel to numeric
countries$PMLevel <- round(as.numeric(countries$PMLevel),2)# Previewing the data
head(countries)## # A tibble: 6 x 2
## Country PMLevel
## <chr> <dbl>
## 1 Estonia 11.1
## 2 Mauritius 11.6
## 3 Australia 13.2
## 4 New Zealand 15
## 5 Ireland 15.2
## 6 Luxembourg 17.5#
countries$iso3c <- as.factor(countrycode(countries$Country, "country.name","iso3c"))url <- "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_0_countries.zip"
folder <- getwd() #set a folder where to download and extract the data
file <- basename(url)
download.file(url, file)
unzip(file, exdir = folder)
#And read it with rgdal library
world <- readOGR(dsn = folder,
layer = "ne_50m_admin_0_countries",
encoding = "UTF-8",
verbose = FALSE)world <- sp::merge(world, countries,
by.x = "ISO_A3",
by.y = "iso3c",
sort = FALSE, duplicateGeoms = TRUE)# Picking a color
pal <- colorNumeric(
palette = "Reds",
domain = countries$PMLevel)# Creating a pop-up for info on a country
world_popup <- paste(world$ADMIN, ",PM10 Level:",world$PMLevel, sep = "")# Adding the map
leaflet(data = world) %>%
addTiles() %>%
setView(0, 0, zoom = 2) %>%
addPolygons(fillColor = ~pal(world$PMLevel),
fillOpacity = 1,
color = "#000000",
weight = 1,
label = ~world_popup) %>%
addLegend("bottomright", pal = pal, values = ~PMLevel, title = "Amount of PM10 Level", opacity =1)