air
Abbigale Mckay
4/30/2020
PM10: Data Retrieval and Processing
# install.packages("rgdal")
# install.packages("ggmap")
# install.packages("leaflet")
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(rgdal)## Loading required package: sp## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
## Linking to sp version: 1.3-2library(countrycode)
library(ggmap)## Loading required package: ggplot2## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.## Please cite ggmap if you use it! See citation("ggmap") for details.library(readxl)
library(leaflet)
library(DT)
library(taucharts)
library(sp)Lets deal with just the PM10 data from sheet 4.
cities<-read_excel("OAP_database.xls", sheet = 4)## New names:
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * `` -> ...6
## * ...cities<-cities[-c(1:2),c(3,2,4)] # some elimination of rows and retention of only 3 columns
names(cities)<-c("City","Country","PMLevel") # Renaming 3 columns
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" ...#worst 100 cities with high levels of PM2.5
# cities$PMLevel<-round(as.numeric(cities$PMLevel),2)
cities$PMLevel<-as.numeric(cities$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: num [1:1099] 42 38 13.5 18.2 13.8 ...cities$PMLevel<-round(cities$PMLevel,2)
cities%>%arrange(-PMLevel)%>%top_n(100)%>%
datatable(cities,rownames=FALSE)## Selecting by PMLevel#which countries contribute the highest number of cities in the top 100 on PM10
tmp<-cities%>%arrange(-PMLevel)%>%top_n(100)%>%
group_by(Country)%>%summarise(number_of_cities=length(Country))## Selecting by PMLeveltmp<-tmp%>%arrange(-number_of_cities)Bar chart of countries with top 100 cities
tmp$Country<-forcats::fct_inorder(tmp$Country)
tauchart(tmp)%>%tau_bar("number_of_cities", "Country",horizontal = "TRUE")%>% tau_legend()%>% tau_tooltip()## Neither color nor size aesthetics have been mapped. Legend plugin will be active but not displayed.#Geographic Map of 100 cities with the highest level of PM2.5
cities100<-cities%>%arrange(-PMLevel)%>%top_n(100)## Selecting by PMLevelcities100$CityCountry<-paste(cities100$City,cities100$Country,sep=", ")
# combining the city and country can help with the geocoding of the city, which appears next.
# locs<-geocode(cities100$CityCountry)
# save(locs,file="locs.RDA")
load("locs.RDA")
cities100$lat<-locs$lat
cities100$lon<-locs$lon## using ggplot2
worldmap <- borders("world", fill="lightgray", colour = "white")
worldmap <- ggplot() + worldmap
worldmap <- worldmap + geom_point(data=cities100,aes(x=lon,y=lat,size=PMLevel),alpha=.4,color="red")
worldmap + theme_void()
The Map Interactive
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)## Assuming "lon" and "lat" are longitude and latitude, respectively#PM 10 countires
countries<-read_excel("OAP_database.xls", sheet = 5)## New names:
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * `` -> ...6
## * ...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.5countries$iso3c=countrycode(countries$Country,"country.name","iso3c")
countries$iso3c<-as.factor(countries$iso3c)#### download the data for the map from naturalearthdata.com
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)
#Colour palette. Check ?colorQuantile and ?RcolorQuantile for more.
pal <- colorNumeric(
palette = "Reds",
domain = countries$PMLevel
)
# Pop up: the information displayed when you click on a country
world_popup <- paste0(world$ADMIN,", PM10 Level:",world$PMLevel,sep = "")
# And 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)