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
I can only imagine how many mouse clicks that would be in a GIS program.
Addendum
- This comment has links to where I acquired the shapefiles (good q by Michael)
- Thanks to this comment there’s now code in a gist that makes a map re-cenetered on the Pacific
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
;-)Can not find the three files. :( the link provided is invalid.
Updated the post. Thx for catching that. This is also the link: http://rud.is/dl/quakefiles.tgz
Pingback: Replicating NatGeo’s “Proper” Earthquake Map in R | Mubashir Qasim
silly q. Where did you get the shapefiles to begin with? Amazing visualization!
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 tozip
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 tocurl
/wget
ordownload.file
that).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
SpatialPointsDataFrame
s 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!
thanks
Pingback: Distilled News | Data Analytics & R
Pingback: Replicating NatGeo’s “Proper” Earthquake Map in R | rud.is
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!!!
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.sorry about the broken code block — here’s a link to a gist https://gist.github.com/achubaty/fe3d6481f152bb14cec4
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.