Bay-SC.1975-2025

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(sf)

## Warning: package 'sf' was built under R version 4.5.2

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(zip)

## 
## Attaching package: 'zip'

## The following objects are masked from 'package:utils':
## 
##     unzip, zip

# Load the data
quakes <- read.csv("Bay-SC.1975-2025.csv")
quakes$time <- as.POSIXct(quakes$time)
quakes$time_numeric <- as.numeric(quakes$time)

# Load faults from KMZ
unzip("qfaults.kmz", exdir = tempdir())
faults <- st_read(file.path(tempdir(), "doc.kml"))

## Multiple layers are present in data source /private/var/folders/l8/gpvbxn2j5dggc88kzrfyy7pr0000gn/T/RtmpRRtkpA/doc.kml, reading layer `Historic Faults'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## automatically selected the first layer in a data source containing more than
## one.

## Reading layer `Historic Faults' from data source 
##   `/private/var/folders/l8/gpvbxn2j5dggc88kzrfyy7pr0000gn/T/RtmpRRtkpA/doc.kml' 
##   using driver `KML'
## Simple feature collection with 14260 features and 2 fields
## Geometry type: MULTILINESTRING
## Dimension:     XYZ
## Bounding box:  xmin: -180 ymin: 18.89354 xmax: 180 ymax: 63.53205
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

# Get quakes boundaries
lon_range <- range(quakes$longitude)
lat_range <- range(quakes$latitude)

# Crop faults to quakes extent
faults <- st_crop(faults, xmin = lon_range[1], xmax = lon_range[2], ymin = lat_range[1], ymax = lat_range[2])

## st_as_s2(): dropping Z and/or M coordinate

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

# Create yearly tick positions from 1975 to 2025
years <- seq(1975, 2025, by = 5)
year_dates <- as.POSIXct(paste0(years, "-01-01"))
year_numerics <- as.numeric(year_dates)

# Start with earthquakes
p <- plot_ly()

# Add fault lines first
for (i in seq_len(nrow(faults))) {
  coords <- st_coordinates(st_geometry(faults[i, ]))
  p <- add_trace(p, x = coords[, 1], y = coords[, 2], z = rep(0, nrow(coords)), 
                 type = 'scatter3d', mode = 'lines',
                 line = list(color = 'red', width = 6), showlegend = (i == 1), 
                 name = 'Faults', hoverinfo = 'skip')
}

# Add earthquake points
p <- add_trace(p, data = quakes, x = ~longitude, y = ~latitude, z = ~depth, 
               type = 'scatter3d', mode = 'markers',
               marker = list(size = ~mag^2, color = ~time_numeric, colorscale = "Rainbow", 
                           opacity = 0.4, showscale = TRUE, 
                           colorbar = list(title = "Year", tickvals = year_numerics, 
                                         ticktext = as.character(years), len = 0.8, 
                                         thickness = 15, tickfont = list(size = 8), 
                                         titlefont = list(size = 10))),
               text = ~paste("Mag:", mag, "<br>Depth:", round(depth, 1), "km<br>Time:", format(time, "%Y-%m-%d")),
               hoverinfo = 'text', name = 'Earthquakes')

p <- layout(p, title = "Bay-SC Earthquakes (1975-Present) with Faults", 
    scene = list(xaxis = list(title = "Longitude"), yaxis = list(title = "Latitude"), zaxis = list(title = "Depth (km)", autorange = "reversed"), aspectratio = list(x = 8, y = 8, z = 6)))

p