Replicating NatGeo’s “Proper” Earthquake Map in R

I saw this post over at NatGeo over the weekend and felt compelled to replicate this:

with ggplot2.

Three shapefiles later and we have it close enough to toss into a post (and I really don’t believe the continent names are necessary).

library(rgdal)
library(ggplot2)
library(ggthemes)
library(ggalt) # devtools::install_github("hrbrmstr/ggalt")

# grab these from http://rud.is/dl/quakefiles.tgz

world <- readOGR("countries.geo.json", "OGRGeoJSON", stringsAsFactors=FALSE)
plates <- readOGR("plates.json", "OGRGeoJSON", stringsAsFactors=FALSE)
quakes <- readOGR("quakes.json", "OGRGeoJSON", stringsAsFactors=FALSE)

world_map <- fortify(world)
plates_map <- fortify(plates)
quakes_dat <- data.frame(quakes)
quakes_dat$trans <- quakes_dat$mag %% 5

gg <- ggplot()
gg <- gg + geom_cartogram(data=world_map, map=world_map,
                          aes(x=long, y=lat, map_id=id),
                          color="white", size=0.15, fill="#d8d8d6")
gg <- gg + geom_cartogram(data=plates_map, map=plates_map,
                          aes(x=long, y=lat, map_id=id),
                          color="black", size=0.1, fill="#00000000", alpha=0)
gg <- gg + geom_point(data=quakes_dat,
                      aes(x=coords.x1, y=coords.x2, size=trans),
                      shape=1, alpha=1/3, color="#d47e5d", fill="#00000000")
gg <- gg + geom_point(data=subset(quakes_dat, mag>7.5),
                      aes(x=coords.x1, y=coords.x2, size=trans),
                      shape=1, alpha=1, color="black", fill="#00000000")
gg <- gg + geom_text(data=subset(quakes_dat, mag>7.5),
                     aes(x=coords.x1, y=coords.x2, label=sprintf("Mag %2.1f", mag)),
                     color="black", size=3, vjust=c(3.9, 3.9, 5), fontface="bold")
gg <- gg + scale_size(name="Magnitude", trans="exp", labels=c(5:8), range=c(1, 20))
gg <- gg + coord_map("mollweide")
gg <- gg + theme_map()
gg <- gg + theme(legend.position=c(0.05, 0.99))
gg <- gg + theme(legend.direction="horizontal")
gg <- gg + theme(legend.key=element_rect(color="#00000000"))
gg

unnamed-chunk-1-1

I can only imagine how many mouse clicks that would be in a GIS program.

Addendum

Cover image from Data-Driven Security
Amazon Author Page

22 Comments Replicating NatGeo’s “Proper” Earthquake Map in R

  1. Tim Salabim

    Hey,
    the link to the files doesn’t seem to work. I would like to use this example as a show and tell for library(mapview). This would be perfect.

    Cheers
    Tim

    Reply
      1. Tim Salabim

        Sweet! Thanks for the quick reply.
        Here’s a quick and dirty approach using mapview

        plates <- readOGR(“plates.json”, “OGRGeoJSON”, stringsAsFactors=FALSE)
        quakes <- readOGR(“quakes.json”, “OGRGeoJSON”, stringsAsFactors=FALSE)
        mag <- exp(quakes@data$mag) / 100

        mapView(plates, fillOpacity = 0) + mapView(quakes, radius = mag)

        Cheers
        Tim

        Reply
        1. hrbrmstr

          it’s almost not fair that only takes 3 additional R expressions (that is, after all the hard work put into mapView ;-)

          Reply
  2. Pingback: Replicating NatGeo’s “Proper” Earthquake Map in R | Mubashir Qasim

    1. hrbrmstr

      Totally not a silly q (I just ran out of time yesterday to add that).

      The world is from Sir Bostock: https://github.com/mbostock/world-atlas

      Plates from ArcGIS: http://www.arcgis.com/home/item.html?id=f155b76c13c84f62864446847f1ae652, just save the lpk file, change extension to zip and unzip it. I converted it to GeoJSON for easier packaging/use.

      Quake data from the USGS: <http://earthquake.usgs.gov/earthquakes/search/> (queries can be saved as a geojson). This specific query was http://earthquake.usgs.gov/fdsnws/event/1/query.geojson?starttime=2015-01-01%2000:00:00&minmagnitude=5&endtime=2015-10-03%2023:59:59&orderby=time-asc (you’ll want to curl/wget or download.file that).

      Reply
  3. James M.

    Super slick post. Do you know of a way to display the map with the Pacific ocean in the center of the map (instead of the Atlantic)? Not much seismicity in the mid-Atlantic. All the good stuff is happening in/around the Pacific.

    Reply
      1. James M.

        Really awesome. Much easier to see the bulk of seismic activity in the re-centered projection. Thanks for sharing the code!

        Reply
  4. Pingback: Distilled News | Data Analytics & R

  5. Pingback: Replicating NatGeo’s “Proper” Earthquake Map in R | rud.is

  6. carloslib

    Sir,

    Excellent post! How do you go by converting the .lpk file into the GeoJSON after you zip it and unzip it?

    Reply
  7. amc

    This is really great, thanks for the post. You inspired me to create a leaflet + shiny version:

    r
    library(leaflet)
    library(shiny)
    leaflet() %>% addTiles() %>%
    addPolygons(data = plates, popup = plates$PLATE,
    color = "black", weight = 1, fillColor = "#00000000") %>%
    addCircleMarkers(lng = quakes_dat$coords.x1, lat = quakes_dat$coords.x2,
    radius = quakes_dat$mag, popup = quakes_dat$title,
    color = "#d47e5d")

    I’m using the default world map from the addTiles() function, and leaflet only does maps in the one projection, so it’s not identical to what you’ve created; however, it’s fun to zoom in an examine the records for individual quakes.

    Reply
  8. Bob

    Thank you so much for this.

    I also tried to understand your “recenter-code”. However, I get this error
    Error in nowrapSpatialPolygons(obj, offset = offset, eps = eps, avoidGEOS = avoidGEOS) :
    unknown coordinate reference system

    when I try
    world <- nowrapRecenter(world)
    .
    How could this be fixed?

    Reply
    1. hrbrmstr

      Take a look at the post and the gist again. I updated the code to work with the new ggplot2 2.x and added some missing library() calls to the gist.

      Reply

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.