knitr::opts_chunk$set(message = FALSE)
library(dplyr)
library(rgdal)
library(countrycode)
library(ggmap)
library(readxl)
library(leaflet)
library(DT)
library(taucharts)
library(sp)

Read and Clean Data

cities=read_excel("OAP_database.xls", sheet=4)
cities = cities[-c(1:2),c(3,2,4)]
names(cities) = c('City','Country','PMLevel')
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" ...

100 Cities with Highest PM2.5 Levels

cities$PMLevel<- round(as.numeric(cities$PMLevel),2)

cities%>%arrange(-PMLevel)%>%top_n(100)%>%datatable(cities,rownames = FALSE)

Countries with the Most Cities in the Top 100 PM10 Pollutors

tmp=cities%>%arrange(-PMLevel)%>%top_n(100)%>%group_by(Country)%>%summarise(number_of_cities=length(Country))%>%arrange(-number_of_cities)

Bar Chart of Countries with Most Top 100 PM10 Cities

tmp$Country=forcats::fct_inorder(tmp$Country)
tauchart(tmp)%>%tau_bar("number_of_cities","Country", horizontal = 'TRUE')%>%tau_legend()%>%tau_tooltip()

Geographic Map of 100 Highest PM2.5 Cities

cities100<-cities%>%arrange(-PMLevel)%>%top_n(100)
cities100$CityCountry=paste(cities100$City,cities100$Country, sep=", ")
locs = geocode(as.character(cities100$CityCountry))

cities100$lat=locs$lat
cities100$lon=locs$lon

worldmap <- borders("world", fill= "lightgray", color = "white")
## Warning: Duplicated aesthetics after name standardisation: colour
worldmap <- ggplot() + worldmap

worldmap <- worldmap + geom_point(data=cities100, aes(x=lon, y=lat, size=PMLevel), alpha=0.4, color='red')
worldmap + theme_void()

Interactive Map

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)%>%addTiles()%>%
  #addProviderTiles('CartoDB.Positron') %>%
  setView(0,0,zoom=2)%>%
  addCircles(stroke=FALSE,fillOpacity = 0.5, color = 'red', radius = ~PMLevel*1000, popup = ~popup)

Interactive Choropleth of Country PM10 Levels

countries=read_excel("OAP_database.xls", sheet=5)
countries<- countries[-c(1:2),c(2:3)]

names(countries) <- c('Country','PMLevel')
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" ...
countries$PMLevel<-round(as.numeric(countries$PMLevel),2)

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)

pal <- colorNumeric(
  palette = 'Reds',
  domain = countries$PMLevel
)

world_popup <- paste(world$admin, ", PM10 Level:", world$PMLevel, sep = "")

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", opacity =1)