Skip navigation

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

19 Comments

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

    • Updated the post. Thx for catching that. This is also the link: http://rud.is/dl/quakefiles.tgz

      Agreed that an interactive view would also be good (I was toying with doing something in the leaflet htmlwidget package, too).

      • 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

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

  2. Can not find the three files. :( the link provided is invalid.

  3. silly q. Where did you get the shapefiles to begin with? Amazing visualization!

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

    • Thx, James. Yep, but I’ll leave re-projecting it to you :-)

      I had to throw together a re-centering function for SpatialPointsDataFrames so rather than clutter the post, here a gist https://gist.github.com/hrbrmstr/7b779df95932752eae1d with the code and an example plot.

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

  5. thanks

  6. Sir,

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

    • You can do something like ogr2ogr -f GeoJSON -tsrs crs:84 plates.json platpol_proj.shp

      • Thanks a lot!!!

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

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

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


3 Trackbacks/Pingbacks

  1. By Replicating NatGeo’s “Proper” Earthquake Map in R | Mubashir Qasim on 05 Oct 2015 at 10:21 am

    […] article was first published on rud.is » R, and kindly contributed to […]

  2. By Distilled News | Data Analytics & R on 05 Oct 2015 at 12:59 pm

    […] Replicating NatGeo’s “Proper” Earthquake Map in R […]

  3. By Replicating NatGeo’s “Proper” Earthquake Map in R | rud.is on 05 Oct 2015 at 1:03 pm

    […] This comment has links to where I acquired the shapefiles (good q by Michael) […]

Leave a Reply

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