eqr

eqr is a thin R wrapper around the USGS FDSN earthquake web service. It returns results as sf objects so they slot directly into spatial workflows, and handles the USGS 20,000-event cap automatically with built-in chunked pagination.

Installation

# install.packages("pak")
pak::pak("gitlab::mhaffner/eqr")

Quick start

library(eqr)

# M2.5+ events globally for the last 30 days (default extent)
quakes <- get_quakes(
  start_time = Sys.time() - 30 * 86400,
  end_time   = Sys.time()
)

quakes
#> Simple feature collection with 1847 features and 26 fields
#> Geometry type: POINT
#> CRS: WGS 84

Functions

`get_quakes()`

Fetches earthquake events and returns an sf data frame. Each row is one event; the geometry column holds a 3-D point (longitude, latitude, depth in km).

get_quakes(
  start_time,
  end_time   = Sys.time(),
  min_mag    = 2.5,
  min_lat    = 24.6,  max_lat = 50,
  min_lng    = -125,  max_lng = -65,
  order_by   = "time",   # "time" | "time-asc" | "magnitude" | "magnitude-asc"
  limit      = NULL,
  auto_chunk = FALSE,
  chunk_days = 30
)

Key columns returned:

Column	Description
`mag`	Magnitude
`place`	Human-readable location description
`time`	Origin time (milliseconds since epoch)
`depth`	Depth in km (extracted from geometry Z coordinate)
`magType`	Magnitude scale (`ml`, `mw`, `mb`, …)
`type`	Event type (usually `"earthquake"`)
`url`	USGS event page
`status`	`"reviewed"` or `"automatic"`

`count_quakes()`

Hits the USGS /count endpoint with the same parameters as get_quakes(), returning only the integer event count. Use this to check whether a query will exceed the 20,000-event cap before fetching data.

n <- count_quakes(
  start_time = "2025-01-01",
  end_time   = "2026-04-15",
  min_mag    = 2.5
)
n
#> [1] 48203

Filtering

By magnitude

major <- get_quakes(
  start_time = "2026-01-01",
  end_time   = "2026-04-15",
  min_mag    = 5.0
)

By region

# Pacific Northwest
pnw <- get_quakes(
  start_time = "2026-01-01",
  end_time   = "2026-04-15",
  min_lat    = 42, max_lat = 49,
  min_lng    = -124, max_lng = -116
)

# Global
global <- get_quakes(
  start_time = Sys.time() - 7 * 86400,
  end_time   = Sys.time(),
  min_lat    = -90, max_lat = 90,
  min_lng    = -180, max_lng = 180
)

Handling large queries

The USGS silently truncates responses at 20,000 events. Always check the count for queries covering long time windows or low magnitude thresholds:

n <- count_quakes(
  start_time = "2025-01-01",
  end_time   = "2026-04-15",
  min_mag    = 2.5
)
n  # if > 20000, use auto_chunk = TRUE

Set auto_chunk = TRUE to have get_quakes() split the window automatically:

quakes <- get_quakes(
  start_time = "2025-01-01",
  end_time   = "2026-04-15",
  min_mag    = 2.5,
  auto_chunk = TRUE,
  chunk_days = 30    # one request per month
)

Internally, auto_chunk calls count_quakes() first, then divides the time window into chunk_days-day segments and combines the results with rbind().

Working with sf output

Because the result is an sf object, the full sf/dplyr ecosystem works immediately:

library(sf)
library(dplyr)

# CRS is always WGS 84
st_crs(quakes)

# Bounding box
st_bbox(quakes)

# Filter with dplyr
shallow <- quakes |> filter(depth < 20)
recent  <- quakes |> filter(as.Date(datetime) == Sys.Date())

# Reproject
quakes_aea <- st_transform(quakes, 9311)

Visualisation

See the vignettes for complete worked examples:

Getting Started — basic usage, filtering, sorting, sf integration
Visualizing Earthquake Data — static and interactive maps with tmap, magnitude distributions, depth vs. magnitude, aftershock sequences
Interactive Maps with mapgl — GPU-accelerated globe maps with data-driven styling, heatmaps, and a spinning global seismicity view

# Quick static map
library(tmap)
tmap_mode("plot")
tm_shape(quakes) +
  tm_bubbles(size = "mag", fill = "depth",
             fill.scale = tm_scale_continuous(values = "magma"))