What you are looking at

R ships with a dataset called quakes: 1,000 seismic events recorded near Fiji since 1964, each with a location, a magnitude, and a depth. Plotted on a map they trace out one of the most active subduction zones on Earth — the boundary where the Pacific Plate dives beneath the Australian Plate.

The map below encodes three variables at once:

Use the layer control in the top right to switch base maps, and to toggle between the individual events and a clustered view. Click any circle for details.

pal <- colorNumeric("plasma", domain = quakes$depth, reverse = TRUE)

popups <- sprintf(
  "<b>Magnitude %.1f</b><br/>Depth: %d km<br/>Recorded by %d stations<br/>%.2f&deg;S, %.2f&deg;E",
  quakes$mag, quakes$depth, quakes$stations, abs(quakes$lat), quakes$long)

# Note: leaflet(quakes) would auto-fit to the data bounds, but these longitudes
# run past 180 degrees, so that auto-fit wraps the globe the wrong way and zooms
# out to the whole world. Passing the data per-layer keeps setView in control.
leaflet() %>%
  addProviderTiles(providers$Esri.WorldImagery,   group = "Satellite") %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas, group = "Minimal") %>%
  addProviderTiles(providers$OpenStreetMap,        group = "Street") %>%

  addCircleMarkers(
    data = quakes, lng = ~long, lat = ~lat,
    radius = ~(mag - 3.5)^2,          # exaggerate so big quakes stand out
    color = ~pal(depth), fillColor = ~pal(depth),
    stroke = TRUE, weight = 1, fillOpacity = 0.6,
    popup = popups, group = "Individual events") %>%

  addMarkers(
    data = quakes, lng = ~long, lat = ~lat, popup = popups,
    clusterOptions = markerClusterOptions(),
    group = "Clustered") %>%

  addLegend("bottomright", pal = pal, values = quakes$depth,
            title = "Depth (km)", opacity = 0.9) %>%

  addLayersControl(
    baseGroups = c("Satellite", "Minimal", "Street"),
    overlayGroups = c("Individual events", "Clustered"),
    options = layersControlOptions(collapsed = FALSE)) %>%

  hideGroup("Clustered") %>%
  setView(lng = 177, lat = -24, zoom = 5) %>%

  # A map that is still below the fold when the page loads has no layout yet, so
  # it renders at world zoom. Re-assert the view once it actually becomes visible.
  htmlwidgets::onRender("function(el, x) {
     var m = this;
     function fix() { m.invalidateSize(); m.setView([-24, 177], 5); }
     if (window.IntersectionObserver) {
       new IntersectionObserver(function(entries) {
         entries.forEach(function(e) { if (e.isIntersecting) fix(); });
       }).observe(el);
     }
     window.addEventListener('load', fix);
     setTimeout(fix, 600);
   }")

The pattern worth noticing

Zoom in on the colours and a structure appears that a plain scatter of dots would hide. The shallow earthquakes (bright yellow, under ~100 km) cluster along the eastern arc near the Tonga trench. Move west and the events get progressively deeper — orange, then purple, down past 600 km.

That gradient is the subducting slab itself. The Pacific Plate meets the trench in the east and descends at an angle to the west, so the further west an earthquake occurs, the deeper down the slab it originated. The map is, in effect, a top-down X-ray of a tectonic plate sliding into the mantle.

par(mar = c(4, 4, 2, 1))
plot(quakes$long, -quakes$depth, pch = 19, cex = 0.5,
     col = pal(quakes$depth),
     xlab = "Longitude (degrees east)", ylab = "Depth (km)",
     main = "Side view: the slab descending from east to west")

Seen in cross-section the slab is unmistakable — a band of earthquakes dipping from near the surface on the right down to 600+ km on the left.

A note on reproducibility

Everything here comes from datasets::quakes, which is built into base R, so this page rebuilds from source with no downloads and no API keys. The map is a Leaflet widget generated by the leaflet R package and embedded directly in the HTML — it is fully interactive offline.

Data: Harvard PRIM-H project, via R’s datasets package. Built 29 September 2026 at 22:37.