library(httr)
library(jsonlite)
library(maps)
library(spatstat.geom)
library(mgcv)

## --- data ---

base <- "https://earthquake.usgs.gov/fdsnws/event/1/query"
q <- list(
  format = "geojson",
  starttime = "2025-07-01",
  endtime   = "2026-02-01",
  minlatitude = 32.4, maxlatitude = 42.1,
  minlongitude = -124.6, maxlongitude = -114.1,
  minmagnitude = 2.5,
  orderby = "time-asc",
  limit = 20000
)
h <- handle(base)
res <- GET(base, query = q, handle = h)
stop_for_status(res)
x <- content(res, as = "text", encoding = "UTF-8")
gj <- fromJSON(x, simplifyVector = FALSE)

feat <- gj$features
df <- data.frame(
  time = as.POSIXct(vapply(feat, \(f) f$properties$time, numeric(1)) / 1000,
                    origin = "1970-01-01", tz = "UTC"),
  mag  = vapply(feat, \(f) f$properties$mag,  numeric(1)),
  lon  = vapply(feat, \(f) f$geometry$coordinates[[1]], numeric(1)),
  lat  = vapply(feat, \(f) f$geometry$coordinates[[2]], numeric(1))
)
df <- df[is.finite(df$lon) & is.finite(df$lat) & !is.na(df$time), ]

## --- spatial plots ---

m_cut <- 4.0
pad <- 2
xlim <- range(df$lon, na.rm = TRUE) + c(-pad, pad)
ylim <- range(df$lat, na.rm = TRUE) + c(-pad, pad)

map("state", fill = FALSE, col = "grey40", lwd = 1, xlim = xlim, ylim = ylim)
points(df$lon, df$lat, pch = 16, col = "grey20",
       cex = pmax(0.4, (df$mag - min(df$mag, na.rm = TRUE) + 0.5) / 2))

## --- nonhomogeneous spatial Poisson via GAM (binned Poisson regression) ---
# observation window
W <- owin(xrange = xlim, yrange = ylim)

# bin the window into a grid
nx <- 20
ny <- 20
xbreak <- seq(W$xrange[1], W$xrange[2], length.out = nx + 1)
ybreak <- seq(W$yrange[1], W$yrange[2], length.out = ny + 1)

# overlay grid used for spatial binning (xbreak, ybreak)
abline(v = xbreak, col = "grey85", lwd = 0.6)
abline(h = ybreak, col = "grey85", lwd = 0.6)

# aggregate within bins
ix <- findInterval(df$lon, xbreak, rightmost.closed = TRUE)
iy <- findInterval(df$lat, ybreak, rightmost.closed = TRUE)
keep <- ix >= 1 & ix <= nx & iy >= 1 & iy <= ny
ix <- ix[keep]; iy <- iy[keep]
y_counts <- as.vector(table(factor(ix, levels = 1:nx),
                            factor(iy, levels = 1:ny)))

# cell midpoints and areas (degrees^2; OK for demo within a small-ish region)
xmid <- (xbreak[-1] + xbreak[-length(xbreak)]) / 2
ymid <- (ybreak[-1] + ybreak[-length(ybreak)]) / 2
grid <- expand.grid(ix = 1:nx, iy = 1:ny)
grid$x <- xmid[grid$ix]
grid$y <- ymid[grid$iy]

dx <- diff(xbreak)[1]
dy <- diff(ybreak)[1]
cell_area <- dx * dy

# arrange in data frame for model input
dat <- data.frame(
  y = y_counts,
  x = grid$x,
  ycoord = grid$y,
  area = cell_area
)

# log E[count] = log(area) + s(x,y)
fit_sp_gam <- gam(y ~ s(x, ycoord, bs = "tp", k = 80) + offset(log(area)),
                  family = poisson(),
                  method = "REML",
                  data = dat)

# plot fitted intensity surface
nx_f <- 200
ny_f <- 200
xg <- seq(W$xrange[1], W$xrange[2], length.out = nx_f)
yg <- seq(W$yrange[1], W$yrange[2], length.out = ny_f)
grid_f <- expand.grid(x = xg, ycoord = yg, area = 1)
lambda_f <- predict(fit_sp_gam, newdata = grid_f, type = "response")
z <- matrix(lambda_f, nrow = nx_f, ncol = ny_f, byrow = FALSE)
zlog <- sqrt(z)
zlim <- quantile(zlog, c(0.01, 0.99), na.rm = TRUE)
image(xg, yg, pmin(pmax(zlog, zlim[1]), zlim[2]),
      col = hcl.colors(128, 
                       palette = "YlOrRd", 
                       rev = TRUE, 
                       alpha = 0.5),
      xlab = "lon", ylab = "lat",
      main = "Fitted spatial intensity (log scale; fine grid)")
map("state", add = TRUE, fill = FALSE, col = "grey40", lwd = 1)
points(df$lon, df$lat, pch = 16, col = "grey20",
       cex = pmax(0.4, (df$mag - min(df$mag, na.rm = TRUE) + 0.5) / 2))