Downloading the data

The data we will be using comes from the US Geological Survey. They have a web service where you can input certain parameters and obtain the data you need. To use the service, click on this link https://earthquake.usgs.gov/earthquakes/search/

In keeping with Dr. Peng’s principle of making a presentation reproducible, the parameters I used for searching for the data are the following:

I also drew a rectangle on the Visayas area of the philippine map but I wasn’t able to log the longitude and latitude values created. I clicked on Search and downloaded the .csv file output of the search. In a future example, i will try using the USGS application programming interfaces (API).

### make sure query.csv file is in your working directory

Oct_2015_Quake <- read.csv("query.csv",               ### read in query.csv in R and assign to R object Oct_2015_Quake
                           stringsAsFactors = FALSE)  ### read factors as characters
Oct_2015_Quake %>% select(time,                       ### select the variables time, latitude, longitude
                          latitude,                   ### mag, and place
                          longitude,
                          mag, place) %>%
        glimpse()                                     ### see data strcuture

## Observations: 83
## Variables: 5
## $ time      <chr> "2013-10-16T23:37:28.370Z", "2013-10-16T22:19:12.410...
## $ latitude  <dbl> 9.7020, 9.9295, 9.7679, 9.8926, 9.8188, 9.5926, 9.48...
## $ longitude <dbl> 123.7936, 123.9820, 123.5686, 124.0694, 123.8910, 12...
## $ mag       <dbl> 5.3, 4.9, 4.8, 4.8, 4.5, 4.9, 4.6, 5.2, 4.9, 4.6, 5....
## $ place     <chr> "6km SW of Bood, Philippines", "2km S of Panaytayon,...

Each row in the dataset contain a recorded earthquake and each variable describes the earthquake in terms of time, location, magnitude and other parameters. Since I am only interested in the three variables i just enumerated I will not bother mentioning the other variables included in the downloaded file.

Tidying data

we see the time variable is read as a string of charactrers. We will transform it to a POSIXct object so we can use it in case we need to denote the time interval between recorded earthquakes. Furthermore, the data uses UTC or Coordinated Universal Time. UTC is the time standard for which the world regulates clocks and time. Philippine Standard Time (PST) is eight hours ahead of UTC. Since R does not recognize PST as a time zone, we will use Asia/Taipei instead.

We will also create another variable concatenating the magnitude of the earthquake and its location.

Oct_2015_Quake <- Oct_2015_Quake %>%                         ### Take the R object we created Oct_2015_Quake
        mutate(date_time = ymd_hms(time,                     ### create new variable from time of class POSIXct
                                   tz = "Asia/Taipei")) %>%  ### Convert from UTC to Asia/Taipei time zone
        mutate(date = date(date_time),                       ### create date variable
               hr = hour(date_time),                         ### create hour variable
               mins = minute(date_time),                     ### create minutes variable
               secs = second(date_time)) %>%                 ### create seconds varaible
        select(date_time,                                    ### selct only the following variables: latitude.
               latitude,                                     ### longitude, mag, place, date, hr, mins, secs
               longitude,
               mag,
               place,
               date,
               hr,
               mins,
               secs)

Oct_2015_Quake$place <- gsub(", Philippines",                ### remove delete ", Philippines" in all rows of 
                             "", Oct_2015_Quake$place)       ### variable Oct_2015$place

Oct_2015_Quake <- Oct_2015_Quake %>%                         ### take R object Oct_2015_Quake
        mutate(mag_loc = paste0(mag,                         ### create new variable mag_loc by combining
                                "Magnitude ",                ### the varaible mag, the strin Magnitude and
                                place))                      ### the variable place

max_mag <- Oct_2015_Quake %>%                                ### take R object Oct_2015_Quake
        dplyr::filter(mag == max(mag))                       ### extract row with highest magnitude earthquake

Oct_2015_Quake %>% select(date_time,                         ### select the following variables and show the
                          latitude,                          ### first six rows
                          longitude,
                          mag,
                          place,
                          mag_loc) %>%
        head()

##             date_time latitude longitude mag                place
## 1 2013-10-17 07:37:28   9.7020  123.7936 5.3       6km SW of Bood
## 2 2013-10-17 06:19:12   9.9295  123.9820 4.9  2km S of Panaytayon
## 3 2013-10-17 02:57:59   9.7679  123.5686 4.8 3km ENE of Dalaguete
## 4 2013-10-17 01:04:31   9.8926  124.0694 4.8   3km SW of Sagbayan
## 5 2013-10-16 23:35:33   9.8188  123.8910 4.5 4km NNW of Antequera
## 6 2013-10-16 23:32:45   9.5926  123.7010 4.9       2km W of Duljo
##                             mag_loc
## 1       5.3Magnitude 6km SW of Bood
## 2  4.9Magnitude 2km S of Panaytayon
## 3 4.8Magnitude 3km ENE of Dalaguete
## 4   4.8Magnitude 3km SW of Sagbayan
## 5 4.5Magnitude 4km NNW of Antequera
## 6       4.9Magnitude 2km W of Duljo

Now, we are ready to plot the data.

Visayas Map

The earthquake occurred in a group of islands called Visayas in the Philippines. The epicenter was located in the island of Bohol 4 Km SE of the town of Sagbayan.

bohol2 <- leaflet() %>%
        setView(lng = 124.1167, 
                lat =  9.8796,
                zoom = 6)
bohol2 %>%
        addTiles() %>%
        addProviderTiles("Esri.NatGeoWorldMap") %>%
        addPopups(lng = max_mag$longitude,
                  lat = max_mag$latitude,
                  popup = "4km SE of Sagbayan, Bohol",
                  options = popupOptions(closeButton = FALSE))

Bohol map

After the initial tremor on tuesday, 82 aftershocks that are greater than 2.5 M have been recorded as of thursday, Oct 17.

icon.pop <- pulseIcons(color = ifelse(Oct_2015_Quake$mag > 7,
                                      'red', ifelse(Oct_2015_Quake$mag < 5,
                                      'palevioletred', 'violetred')),
                       heartbeat = ifelse(Oct_2015_Quake$mag > 7,
                                          0.4, ifelse(Oct_2015_Quake$mag < 5,
                                      1, 0.8)),
                       iconSize = ifelse(Oct_2015_Quake$mag > 7,
                                          12, ifelse(Oct_2015_Quake$mag < 5,
                                      4, 8)))
bohol3 <- leaflet() %>%
        setView(lng = 124.1167, 
                lat =  9.8796,
                zoom = 9)

bohol3 %>%
        addTiles() %>%
        addProviderTiles("OpenStreetMap.Mapnik") %>%
        addPulseMarkers(
    lng = Oct_2015_Quake$longitude, lat = Oct_2015_Quake$latitude,
    label= Oct_2015_Quake$place,
    icon = icon.pop)

Reading in Typhoon Data

Data for typhoon haiyan came from the National Oceanic and Atmospheric Administration (NOAA). The International Best Track Archive for Climate Stewardship (IBTrACS) is a global set of historical tropical cyclones. The data is derived from contributions from various world wide agencies tasked to monitor climate and weather. Description of the data and instructions for downloading the data are available at https://www.ncdc.noaa.gov/ibtracs/index.php?name=ibtracs-data-access

The storm data for the year 2013 was downloaded from this ftp site. The variable names contain prefixes to denote the institute/agency that contributed the data. For our purposes we will use the data from the Joint typhoon Warning System (JWTC).

ibtracks2013 <- read.csv("Year.2013.ibtracs_all.v03r10.csv", ### take the downloaded file from NOAA and read
                         skip = 1,                           ### content but skip firsst row
                         stringsAsFactors = FALSE,           ### read factors as strings or characters
                         na.strings = "-999.0")              ### specify code for missing values

var_names <- c("Name", "ISO_time")                           ### create var names containing Name and ISO_time

haiyan <- ibtracks2013 %>%                                   ### take the r object ibtracks2013 we created
        filter(Name == "HAIYAN") %>%                         ### filter typhoon haiyan
        select(var_names, starts_with("jtwc"))               ### select varaibles Name, ISO_time, and all variables beginning with the string jtwc

jtwc <- haiyan %>%                                           ### take r object haiyan and
        select(starts_with("jtwc"))                          ### select only variables beginning with jtwc

non_jtwc <- haiyan %>%                                       ### take r object haiyan and
        select(1:4)                                          ### select only variables Name, ISO_time, jtwc_sh_lat

non_jtwc$ISO_time <- ymd_hms(non_jtwc$ISO_time)            ### convert ISO_time to POSIXct object         

jtwc_df <- as.data.frame(rbind(sapply(jtwc, as.numeric)))    ### convert all variables in jtwc_df to numeric and row bind

var_empty <- sapply(jtwc_df, function(x) sum(is.na(x)))      ### check which variables are empty/ all missing data

index <- var_empty == 38                                     ### create index which variables contain only NAs 

plot_df <-cbind(non_jtwc, jtwc_df[, !index])                 ### use index to remove varaibles that contain only NAs 

plot_df <- plot_df %>%                                       ### remove last row
        select(-starts_with("jtwc_sh")) %>% 
        slice(1:37)

plot_df                                                      ### show first 10 rows of data

## # A tibble: 37 x 22
##    Name   ISO_time            jtwc_wp_lat jtwc_wp_lon jtwc_wp_wind
##    <chr>  <dttm>                    <dbl>       <dbl>        <dbl>
##  1 HAIYAN 2013-11-02 06:00:00         6.6        162.           15
##  2 HAIYAN 2013-11-02 12:00:00         6.2        160.           15
##  3 HAIYAN 2013-11-02 18:00:00         6.1        159.           15
##  4 HAIYAN 2013-11-03 00:00:00         6.2        157.           20
##  5 HAIYAN 2013-11-03 06:00:00         6.3        156.           25
##  6 HAIYAN 2013-11-03 12:00:00         6.6        154.           30
##  7 HAIYAN 2013-11-03 18:00:00         6.3        153.           30
##  8 HAIYAN 2013-11-04 00:00:00         6.1        152.           35
##  9 HAIYAN 2013-11-04 06:00:00         6.1        150.           40
## 10 HAIYAN 2013-11-04 12:00:00         6.3        149.           45
## # ... with 27 more rows, and 17 more variables: jtwc_wp_pres <dbl>,
## #   jtwc_.._rmw <dbl>, jtwc_.._poci <dbl>, jtwc_.._roci <dbl>,
## #   jtwc_.._eye <dbl>, jtwc_.._wrad34_rad1 <dbl>,
## #   jtwc_.._wrad34_rad2 <dbl>, jtwc_.._wrad34_rad3 <dbl>,
## #   jtwc_.._wrad34_rad4 <dbl>, jtwc_.._wrad50_rad1 <dbl>,
## #   jtwc_.._wrad50_rad2 <dbl>, jtwc_.._wrad50_rad3 <dbl>,
## #   jtwc_.._wrad50_rad4 <dbl>, jtwc_.._wrad64_rad1 <dbl>,
## #   jtwc_.._wrad64_rad2 <dbl>, jtwc_.._wrad64_rad3 <dbl>,
## #   jtwc_.._wrad64_rad4 <dbl>

Description of the Variables

Below is a brief description of the data from the Joint Typhoon Warning Center (JTWC).

Tidying the data

We will need to create shorter and more decriptive names to our variables and organize them in a way that will make plotting easier.

names(plot_df) <- c("name", "date_time", "latitude",       ### change variable names to more meaningful names
                    "longitude", "wind", "pressure",
                    "max_wind_rad", "poci", "roci",
                    "eye_diameter", "34_ne", "34_se",
                    "34_sw", "34_nw", "50_ne",
                    "50_se", "50_sw", "50_nw",
                    "64_ne", "64_se", "64_sw", "64_nw")

plot_df <- plot_df %>% select(name,                       ### select variables name, date_time, latitude, longitude
                              date_time,                  ### and variables 11 to 22
                              latitude,
                              longitude,
                              11:22) %>%
        gather(key = wind_spd_dir,                        ### gather all varaible to column wind_spd_dir except variables
               value = value,                             ### name, date_time, latitude, longitude
               -name,
               -date_time,
               -latitude,
               -longitude,
               na.rm = TRUE) %>%                          ### don't include NAs or missing values
        separate(col = wind_spd_dir,                      ### separate wind_spd_dir to wind_spd and dir
                 c("wind_spd", "dir")) %>%
        spread(dir, value) %>%
        arrange(date_time,                                ### arrange rows according to date_time, latitude, 
                latitude,                                 ### and longitude
                longitude)

plot_df$wind_spd <- as.factor(plot_df$wind_spd)           ### convert variable wind_spd to factor

nov4AM6 <-plot_df %>% slice(25:27)                        ### create dataframe with only row 25:27
nov5AM6 <-plot_df %>% slice(37:39)
nov6AM6 <-plot_df %>% slice(49:51)
nov7AM6 <-plot_df %>% slice(61:63)
nov8AM0 <-plot_df %>% slice(79:81)
nov8bohol <-plot_df %>% slice(70:72)
nov9AM18 <-plot_df %>% slice(91:93)
nov10AM12 <-plot_df %>% slice(100:102)
nov11AM6 <-plot_df %>% slice(109:111)

nov8bohol

## # A tibble: 3 x 9
##   name   date_time           latitude longitude wind_spd    ne    nw    se
##   <chr>  <dttm>                 <dbl>     <dbl> <fct>    <dbl> <dbl> <dbl>
## 1 HAIYAN 2013-11-08 00:00:00       11      125. 34         130   130   115
## 2 HAIYAN 2013-11-08 00:00:00       11      125. 50          65    70    60
## 3 HAIYAN 2013-11-08 00:00:00       11      125. 64          50    50    45
## # ... with 1 more variable: sw <dbl>

For this part, I will be plotting the data using the geom from the Building Data Visulaization Course I just finished. The data is organized into a dataframe consisting of three rows. Each row contains the wind speed (35, 50 and 64 knots) in which the data is reported for a unique latitude and longitude location. Each row also contains the radius of the wind speed in 4 directions (northeast, southeast, southwest and northwest).

### R code to create geom_hurricane

GeomHurricane <- ggplot2::ggproto(
        `_class` = "GeomHurricane",          ### name of new geom
        `_inherit` = Geom,                   ### class of R object
        required_aes = c("x",                ### argument for lingitude
                         "y",                ### argument for latitude
                         "r_ne",             ### argument for NE wind radius
                         "r_se",             ### argument for SE wind radius
                         "r_sw",             ### argument for SW wind radius
                         "r_nw"),            ### argument for NW wind radius
        default_aes = aes(colour = "NA",     ### default value for colour argument
                          fill = "NA",       ### default value for fill argument
                          linetype = 1,      ### default value for linetype argument
                          alpha = 0.8,       ### default value for alpha argument
                          scale_radii = 1),  ### default value for scale_radii argument
        draw_key = draw_key_polygon,         ### setting draw_key_polygon to draw legend
        draw_group = function(data,                 ### selecting draw_group over draw_panel
                              panel_scales,         ### elements are grouped according to wind_speed
                              coord){              ### data describing coordinates
                center_pt <- c(data[1, ]$x,         ### coordinates for longitude
                               data[1, ]$y)         ### coordinates for latitude
                color <- data[1, ]$colour           ### selecting colour from first row
                fill <- data[1, ]$fill              ### selecting fill from first row
                alpha <- data[1, ]$alpha            ### selecting alpha from first row
                scale_radii = data[1, ]$scale_radii ### selecting scale_radii from first row
                ne_qtr_circ <- geosphere::destPoint(p = center_pt, ### setting center of arc
                                                    b = 0:90,      ### setting points for arc
                                                    d = data[1, ]$r_ne * 1852 * scale_radii) ### setting radius
                ### converted nautical miles to meters by multiplying aby 1852 and to map scale by scale_radii
                
                ne_df <- data.frame(x = c(ne_qtr_circ[, "lon"],
                                          center_pt[1]),
                                    y = c(ne_qtr_circ[, "lat"],
                                          center_pt[2])
                )
                se_qtr_circ <- geosphere::destPoint(p = center_pt, ### setting center of arc
                                                    b = 90:180,      ### setting points for arc
                                                    d = data[1, ]$r_se * 1852 * scale_radii) ### setting radius
                ### converted nautical miles to meters by multiplying aby 1852 and to map scale by scale_radii
                
                se_df <- data.frame(x = c(se_qtr_circ[, "lon"],
                                          center_pt[1]),
                                    y = c(se_qtr_circ[, "lat"],
                                          center_pt[2])
                )
                
                nw_qtr_circ <- geosphere::destPoint(p = center_pt, ### setting center of arc
                                                    b = 270:360,      ### setting points for arc
                                                    d = data[1, ]$r_nw * 1852 * scale_radii) ### setting radius
                ### converted nautical miles to meters by multiplying aby 1852 and to map scale by scale_radii
                
                nw_df <- data.frame(x = c(nw_qtr_circ[, "lon"],
                                          center_pt[1]),
                                    y = c(nw_qtr_circ[, "lat"],
                                          center_pt[2])
                )
                
                sw_qtr_circ <- geosphere::destPoint(p = center_pt, ### setting center of arc
                                                    b = 180:270,      ### setting points for arc
                                                    d = data[1, ]$r_sw * 1852 * scale_radii) ### setting radius
                ### converted nautical miles to meters by multiplying aby 1852 and to map scale by scale_radii
                
                sw_df <- data.frame(x = c(sw_qtr_circ[, "lon"],
                                          center_pt[1]),
                                    y = c(sw_qtr_circ[, "lat"],
                                          center_pt[2])
                )
                
                all_df <- rbind(ne_df, se_df, nw_df, sw_df)
                coords <- coord$transform(all_df,
                                          panel_scales)
                
                grid::polygonGrob(x = coords$x,
                                  y = coords$y,
                                  gp = grid::gpar(col = color,
                                                  fill = fill,
                                                  alpha =alpha))
        }
                          
)

geom_hurricane <- function(mapping = NULL,
                           data = NULL, 
                           stat = "identity",
                           position = "identity",
                           na.rm = FALSE,
                           show.legend = NA,
                           inherit.aes = TRUE,
                           ...) {
        ggplot2::layer(
                geom = GeomHurricane,
                mapping = mapping,
                data = data,
                stat = stat,
                position = position,
                show.legend = show.legend,
                inherit.aes = inherit.aes,
                params = list(na.rm = na.rm,
                              ...)
        )
}

Hurricane Haiyan over Bohol, Samar, and Leyte

Barely 3 weeks after the 7.1 M earthquake struck at Sagabayan, Bohol, Super Typhoon Haiyan made landfall at the town of Guiuan, Eastern Samar. The red, orange, and yellow area indicating wind speed of 34, 50, and 64 knots winds (63, 93, 118 km/hr).

leyte_casualties <- read.csv("leyte_casualties.csv")
phil_map <- readRDS("PHL_adm2.rds")                           ###read shape files for Philippine Map
bohol_map <- phil_map[phil_map$NAME_1 == "Bohol", ]           ### filter only map for Bohol
bohol_mapdf <- fortify(bohol_map)                             ### use fortify function so you can use map in ggplot

sagbayan_map <- phil_map[phil_map$NAME_2 == "Sagbayan", ]     ### filter map only for sagbayan
sagbayan_mapdf <- fortify(sagbayan_map)                       ### use fortify function so you can use map in ggplot

Esamar_map <- phil_map[phil_map$NAME_1 == "Eastern Samar", ]
Esamar_mapdf <- fortify(Esamar_map)                        


guiuan_map <- phil_map[phil_map$NAME_2 == "Guiuan", ]
guiuan_mapdf <- fortify(guiuan_map)

dist <- mapdist("Sagbayan, Bohol", "Guiuan, Samar")

guiuan <- geocode("Guiuan, Samar")

towns <- data.frame(list(places = c("Guiuan",     ### create dataframe for labels for Guiuan and sagbayan
                                    "Sagbayan",
                                    "Tacloban\nCity"),
                         lon = c(125.9245,
                                 123.9975,
                                 124.7644),
                         lat = c(10.53525,
                                 10.09920,
                                 11.11317)))
prov <- data.frame(list(places = c("East Samar",  ### create dataframe for labels for East Samar and Bohol
                                   "Bohol"),
                        lon = c(125.9950,
                                124.1435),
                        lat = c(11.758273,
                                9.449991)))


leyte_map <- phil_map[phil_map$NAME_1 == "Leyte", ]           ### filter only map for Leyte
leyte_mapdf <- fortify(leyte_map, region = "ID_1")                             ### use fortify function so you can use map in ggplot

tacloban_map <- phil_map[phil_map$NAME_2 == "Tacloban City", ]     ### filter map only for Tacloban
tacloban_mapdf <- fortify(tacloban_map, region = "ID_2")                       ### use fortify function so you can use map in ggplot


leyte_towns_merged <- geo_join(leyte_map, leyte_casualties, "NAME_2", "town")


map_plot <- get_map(c(lon = 124.7,                  ### get map from ggogle map
                      lat = 11),                    ### specify arguments for longitude, latitude,zoom
                    zoom = 8,                       ### and maptype
                    maptype = "toner-background") 

covered <- map_plot %>%                             ### use map we created earlier
  ggmap(extent = "device") +                        ### specify border of map
  geom_hurricane(data = nov8bohol,                  ### specify dataframe to use
                 aes(x = longitude,                 ### specify argument for longitude
                     y = latitude,                  ### specify argument for latitude
                     r_ne = ne,                     ### specify northeast radius
                     r_se = se,                     ### specify southeat radius
                     r_nw = nw,                     ### specify northwest radius
                     r_sw = sw,                     ### specify southwest radius
                     color = wind_spd,              ### specify argument for color
                     fill = wind_spd),              ### specify argument for wind_spd
                 scale_radii = 1) +                 ### specify what percent of radii will be shown
  scale_color_manual(name = "Wind speed (kts)",     ### create legend
                     values = c("red",
                                "orange",
                                "yellow")) + 
  scale_fill_manual(name = "Wind speed (kts)", 
                    values = c("red",
                               "orange",
                               "yellow")) +
        geom_polygon(data = bohol_mapdf,            ### create map of bohol
                     aes(x = long,
                         y = lat,
                         group = group)) +
        geom_polygon(data = sagbayan_mapdf,         ### create map of sagbayan and color it green
                     aes(x =long,
                         y = lat,
                         group = group),
                     fill = "green") +
        geom_point(data = max_mag,                  ### mark epicenter of 7.1 M earthquake
                   aes(x= longitude,
                       y = latitude),
                   color = "blue", size = 2) +
        geom_polygon(data = Esamar_mapdf,           ### create map of eastern samar
                     aes(x= long,
                         y = lat,
                         group = group)) +
        geom_polygon(data = guiuan_mapdf,           ### create map of guiuan and color it green
                     aes(x = long,
                         y = lat,
                         group = group),
                     fill = "green") +
        geom_polygon(data = leyte_mapdf,           ### create map of eastern samar
                     aes(x= long,
                         y = lat,
                         group = group)) +
        geom_polygon(data = tacloban_mapdf,           ### create map of guiuan and color it green
                     aes(x = long,
                         y = lat,
                         group = group),
                     fill = "green") +
        geom_text(data = prov,                      ### create labels for Eastern samar and Bohol
                  aes(x = lon,
                      y = lat,
                      label = places),
                  color = "gray7") +
        geom_text(data = towns,                     ### create labels for guiuan and sagbayan
                  aes(x = lon,
                      y = lat,
                      label = places),
                  color = "green") +
        guides(col=guide_legend(title="Wind\nspeed\n(kts)"))

covered