1) Overview

 

1-1) Motivation for this project

For this project, I was looking for a dataset which I could visualize it on a map. And then I came across the news of an earthquake occurred near Fiji on July 26th. Earthquake seemed like a good dataset to show on map.

 

1-2) Challenges

Fiji sits on the Pacific ‘Ring of Fire’ meaning that earthquakes happen frequently in this region. While having sufficient data is a plus usually, I was concerned whether there will be too much data. In that case, it might be difficult to visualize the graph effectively.

 

1-3) Proposed Design

First, we will show the world map so people who are unfamiliar with the location of Fiji could get the general idea of its location. The map will also show Fiji’s population

Then we will show earthquakes of that region. The data comprises data from 1960 of magnitude higher than 4. This will be a step by step process where users could see thoroughly how it was built.

Finally, we will make some basic data analysis with the data.

 

 

2) Data Exploration

 

2-1) Importing necessary libraries

In this project we will be needing the following libraries

library(shiny)
library(shinydashboard)
library(leaflet)
library(leaflet.extras)
library(rgdal)
library(sp)
library(raster)

 

2-2) Data wrangling

We will be using one of R datasets package called ‘quakes’. All the R datasets can be found here(https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/00Index.html#Q).

 

2-2-1) head(quakes)

Luckily we find that the data consists of geographical information specifically ‘lat’ and ‘long’ which stand for latitude and longitude.

head(quakes)
##      lat   long depth mag stations
## 1 -20.42 181.62   562 4.8       41
## 2 -20.62 181.03   650 4.2       15
## 3 -26.00 184.10    42 5.4       43
## 4 -17.97 181.66   626 4.1       19
## 5 -20.42 181.96   649 4.0       11
## 6 -19.68 184.31   195 4.0       12

 

2-2-2) summary(quakes)

Let’s check the summary to get the idea of the dataset we’re dealing with. We see that ‘mag’ which probably stands for ‘magnitude’. We check that the ‘min’ value of ‘mag’ is 4.

summary(quakes)
##       lat              long           depth            mag      
##  Min.   :-38.59   Min.   :165.7   Min.   : 40.0   Min.   :4.00  
##  1st Qu.:-23.47   1st Qu.:179.6   1st Qu.: 99.0   1st Qu.:4.30  
##  Median :-20.30   Median :181.4   Median :247.0   Median :4.60  
##  Mean   :-20.64   Mean   :179.5   Mean   :311.4   Mean   :4.62  
##  3rd Qu.:-17.64   3rd Qu.:183.2   3rd Qu.:543.0   3rd Qu.:4.90  
##  Max.   :-10.72   Max.   :188.1   Max.   :680.0   Max.   :6.40  
##     stations     
##  Min.   : 10.00  
##  1st Qu.: 18.00  
##  Median : 27.00  
##  Mean   : 33.42  
##  3rd Qu.: 42.00  
##  Max.   :132.00

 

 

3) World Map

Using world map, we can see that Fiji is situated in Oceania. North to New Zealand and East to Papua New Guinea. It’s a small country with only 830,000 people.

download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip" , destfile="world_shape_file.zip")
unzip("world_shape_file.zip")
world_spdf = readOGR(dsn=getwd(), layer="TM_WORLD_BORDERS_SIMPL-0.3")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\User\Desktop\SMU_S2\2.Visualization_Analytics\Assignment\Assignment5\Earthquake", layer: "TM_WORLD_BORDERS_SIMPL-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings:  POP2005
world_spdf$POP2005 = as.numeric(as.character(world_spdf$POP2005)) / 1000000 %>% round(2)


#Choropleth Map
bins=c(0,10,20,50,100,500,Inf)
pal = colorBin(palette = "YlOrBr", domain=world_spdf$POP2005, na.color = "transparent", bins=bins)

customLabel = paste("Country: ", world_spdf$NAME, "<br/>", "Population: ", round(world_spdf$POP2005, 2), sep = "") %>% 
  lapply(htmltools::HTML)

leaflet(world_spdf) %>%
  
  addProviderTiles(providers$OpenStreetMap, options = tileOptions(minZoom=2, maxZoom=8)) %>%
  
  addPolygons(fillColor = ~pal(POP2005),
              fillOpacity = 0.9,
              stroke = TRUE,
              color="white",
              highlight=highlightOptions(
                weight = 5,
                fillOpacity = 0.3
              ),
              label = customLabel,
              weight = 0.3,
              smoothFactor = 0.2) %>%
  
  addLegend(
    pal=pal,
    values = ~POP2005,
    position = "bottomright",
    title ="World Population (Millions)"
  ) %>%
  
  setView(178.0650, -17.7134, zoom = 3.5)

 

 

4) Plot Earthquake

We’ll be explaining step by step of how to plot the earthquake on a map

 

4-1) Show a map

We are going to use OpenStreetMap as our base map

leaflet(data = quakes) %>% 
      addTiles(group = "OpenStreetMap")

 

4-2) Add magnitude

The circles will indicate the magnitude. Size of the circles will show how big the earthquake was.

pal <- colorNumeric("OrRd", quakes$mag)

leaflet(data = quakes) %>% 
  addTiles(group = "OpenStreetMap") %>%
  
  addCircles(radius = ~10^quakes$mag/10, weight=1, color=~pal(quakes$mag), fillColor = ~pal(quakes$mag),
             fillOpacity = 0.7, popup = as.character(quakes$mag), label = ~as.character(quakes$mag), 
             group = "Points")
## Assuming "long" and "lat" are longitude and latitude, respectively

 

4-3) Add a legend

Without a legend, it’s difficult to see how big the earthquake was. So at bottom right, we are going to create a legend of how a smaller and lighter color represent smaller magnitude whereas bigger and darker color represent bigger earthquakes.

pal <- colorNumeric("OrRd", quakes$mag)

leaflet(data = quakes) %>% 
  addTiles(group = "OpenStreetMap") %>%
  
  addCircles(radius = ~10^quakes$mag/10, weight=1, color=~pal(quakes$mag), fillColor = ~pal(quakes$mag),
             fillOpacity = 0.7, popup = as.character(quakes$mag), label = ~as.character(quakes$mag), 
             group = "Points") %>%
  
  
  addLegend(
    position = "bottomright",
    pal = pal,
    values = ~mag, group = "Points",
    title = "Earthquake Magnitude"
  )
## Assuming "long" and "lat" are longitude and latitude, respectively

 

4-4) Add layers control

Although ‘openStreetMap’ serves its purpose, we will include more maps to suit tastes of various people. Users could choose the map that’s most appealing for oneself.

pal <- colorNumeric("OrRd", quakes$mag)

leaflet(data = quakes) %>% 
  
  addTiles(group = "OpenStreetMap") %>%
  
  addProviderTiles(providers$Esri.WorldStreetMap, options = tileOptions(minZoom = 0, maxZoom = 13), group = "Esri.WorldStreetMap") %>%
    
  addProviderTiles(providers$Esri.WorldImagery, options = tileOptions(minZoom = 0, maxZoom = 13), group = "Esri.WorldImagery") %>%
  
  addCircles(radius = ~10^quakes$mag/10, weight=1, color=~pal(quakes$mag), fillColor = ~pal(quakes$mag),
               fillOpacity = 0.7, popup = as.character(quakes$mag), label = ~as.character(quakes$mag), 
               group = "Points") %>%
  
  
  addLegend(
    position = "bottomright",
    pal = pal,
    values = ~mag, group = "Points",
    title = "Earthquake Magnitude"
  ) %>%


  addLayersControl(
    baseGroups = c("OpenStreetMap", "Esri.WorldStreetMap", "Esri.WorldImagery"),
    options = layersControlOptions(collapsed = TRUE)
  )
## Assuming "long" and "lat" are longitude and latitude, respectively

 

4-5) Show number of earthquakes

At the moment, it is difficult to visualize how many earthquakes there were. We will cluster the dots with a number on top of it to get the general idea of how many earthquakes occurred at which area. Since too much data might overwhelm users, we will put this into layer so users could choose whether to see it or not.

pal <- colorNumeric("OrRd", quakes$mag)
eliminate <- quakes[!is.na(quakes$long)&!is.na(quakes$lat),]

leaflet(data = quakes) %>% 
  
  addTiles(group = "OpenStreetMap") %>%
  
  addProviderTiles(providers$Esri.WorldStreetMap, options = tileOptions(minZoom = 0, maxZoom = 13), group = "Esri.WorldStreetMap") %>%
    
  addProviderTiles(providers$Esri.WorldImagery, options = tileOptions(minZoom = 0, maxZoom = 13), group = "Esri.WorldImagery") %>%
  
  addCircles(radius = ~10^quakes$mag/10, weight=1, color=~pal(quakes$mag), fillColor = ~pal(quakes$mag),
               fillOpacity = 0.7, popup = as.character(quakes$mag), label = ~as.character(quakes$mag), 
               group = "Points") %>%
  
  addMarkers(
    popup=as.character(quakes$mag), 
    label=as.character(quakes$mag),
    clusterOptions = markerClusterOptions(),
    group = "Groups"
  ) %>%
  
  
  addLegend(
    position = "bottomright",
    pal = pal,
    values = ~mag, group = "Points",
    title = "Earthquake Magnitude"
  ) %>%


  addLayersControl(
    baseGroups = c("OpenStreetMap", "Esri.WorldStreetMap", "Esri.WorldImagery"),
    overlayGroups = c("Points", "Groups"),
    options = layersControlOptions(collapsed = TRUE)
  )
## Assuming "long" and "lat" are longitude and latitude, respectively
## Assuming "long" and "lat" are longitude and latitude, respectively

 

4-6) Heatmap

We could make an assumption that a heatmap will correspond to the number earthquakes. For example, the more the occurrences, the more ‘red’ it will be represented in the heatmap.

It could give us an idea about where the ‘Ring of Fire’ is.

On R shiny, heatmap will be included separately in a separate sidebar.

# Heatmap
pal <- colorNumeric("RdYlBu", quakes$mag)

leaflet(data = quakes) %>% 
  
  addProviderTiles(providers$Esri.WorldImagery, options = tileOptions(minZoom = 0, maxZoom = 13)) %>%

  addHeatmap(
    ~long, ~lat, intensity = quakes$mag, blur = 20, max = 0.05, radius = 15
  )

 

5) Further analysis

Remember during the data wrangling stage, we were able to understand ‘lat’, ‘long’, ‘depth’, and ‘mag’ but there was one more column called ‘stations’ which we haven’t touched upon yet. It turns out that ‘stations’ concerns the number of stations reporting this earthquake.

colnames(quakes)
## [1] "lat"      "long"     "depth"    "mag"      "stations"

We could think of a simple linear regression. If a high magnitude earthquake occurred, more stations would have reported it to the world. Let’s find out if our thought was correct.

We are going to plot a diagram with X: magnitude and Y: # of stations reported.

When we check the data. We confirmed that our thought was correct.

attach(quakes)
Quake.mod <- lm(stations ~ mag)
Quake.mod
## 
## Call:
## lm(formula = stations ~ mag)
## 
## Coefficients:
## (Intercept)          mag  
##     -180.42        46.28
plot(jitter(mag, amount = 0.05), stations, 
     pch = 20,
     ylab = "# of Stations Reporting",
     xlab = "Magnitude",
     main = "Fiji Earthquakes Magnitude and Reporting",
     col = rgb(0.1, 0.2, 0.8, 0.3))

abline(-180.42, 46.28, col="red", lwd = 2)

6) Final thoughts

R shiny was very challenging.

Through this project, we were able to see that Fiji is indeed along the ‘Ring of Fire’ with many earthquakes.

We found out that ,

Thank you!