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)
 
# 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_map(data=world_map, map=world_map,
                    aes(x=long, y=lat, map_id=id),
                    color="white", size=0.15, fill="#d8d8d6")
gg <- gg + geom_map(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

Buy on AmazonDDS Blog
DDS PodcastAmazon Author Page

20 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. 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

Leave a Reply