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(ggalt) # devtools::install_github("hrbrmstr/ggalt")

# grab these from

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"))


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


Buy on AmazonDDS Blog
DDS PodcastAmazon 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

      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

        1. hrbrmstr

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

  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:

      Plates from ArcGIS:, 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: <> (queries can be saved as a geojson). This specific query was (you’ll want to curl/wget or download.file that).

  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.

      1. James M.

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

  4. Pingback: Distilled News | Data Analytics & R

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

  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.

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

    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.


Leave a Reply