gmsp processes ground-motion time histories, computes
response spectra and intensity measures, selects and matches records,
and composes modal decomposition and reconstruction through
xemd.
Time-history tables use OCID, ID,
t and s for component, quantity, time in
seconds and sample value. The by argument names the
metadata columns that identify a record; character()
declares one record without metadata. The synthetic record below has two
horizontal acceleration components, H1 and H2,
sampled every 0.01 s and identified by RecordID.
library(gmsp)
library(data.table)
Time <- seq(0, 20, by = 0.01)
Envelope <- sin(pi * Time / max(Time))^2
TSL <- rbindlist(list(
data.table(RecordID = "R1", OCID = "H1", ID = "AT", t = Time,
s = 1200 * Envelope * sin(2 * pi * 1.5 * Time)),
data.table(RecordID = "R1", OCID = "H2", ID = "AT", t = Time,
s = 900 * Envelope * sin(2 * pi * 2.5 * Time + 0.4))
))
head(TSL, 3)
#> RecordID OCID ID t s
#> <char> <char> <char> <num> <num>
#> 1: R1 H1 AT 0.00 0.0000000000
#> 2: R1 H1 AT 0.01 0.0002786433
#> 3: R1 H1 AT 0.02 0.0022192480buildTSL() processes one measured quantity into
acceleration (AT), velocity (VT) and
displacement (DT) on a common time grid per record. The
frequency of interest fmax, the STFT window length and
overlap, the edge-taper thresholds astop0 and
apass0 and the low-pass magnitudes astopLP and
apassLP are required arguments. The control values in this
guide illustrate the calls. Raw input is converted from
unitsSource, regularized and tapered, and the amplitudes
are returned in unitsTarget/s^2, unitsTarget/s
and unitsTarget.
Histories <- buildTSL(
TSL, by = "RecordID", unitsSource = "mm", unitsTarget = "mm", fmax = 5,
windowLength = 128, overlap = 75,
astop0 = 1e-3, apass0 = 1e-3,
astopLP = 1e-3, apassLP = 0.98,
includeWarnings = FALSE
)
Histories[, .(samples = .N, start = min(t), peak = max(abs(s))), by = .(OCID, ID)]
#> OCID ID samples start peak
#> <char> <char> <int> <num> <num>
#> 1: H1 AT 798 0 1199.068679
#> 2: H1 VT 798 0 130.016874
#> 3: H1 DT 798 0 14.212090
#> 4: H2 AT 798 0 899.163994
#> 5: H2 VT 798 0 57.569217
#> 6: H2 DT 798 0 3.686879With kNyq = NULL, the default, buildTSL()
selects the sampling frequency; the processed grid of this record has a
step of 0.025 s. Tapering, resampling and the removal of means affect
the processed histories, and the returned quantities are not exact
algebraic inverses of the input.
computeSpectraTSL() computes PSA, PSV and SD from
AT, VT and DT as separate
excitations. Each ordinate retains the units of its source domain, and
the zero period carries the peak absolute value of the regularized
source signal. includeD50 adds D50, the median
over nTheta angles in [0, pi) of the rotated horizontal
component, which requires H1 and H2.
PSL <- computeSpectraTSL(Histories, by = "RecordID", xi = 0.05,
Tn = c(0.1, 0.2, 0.5, 1, 2),
includeD50 = TRUE, nTheta = 8)
PSL[ID == "PSA" & OCID == "D50"]
#> RecordID OCID Tn ID S
#> <char> <char> <num> <char> <num>
#> 1: R1 D50 0.0 PSA 1344.4011
#> 2: R1 D50 0.1 PSA 1389.0485
#> 3: R1 D50 0.2 PSA 1603.9814
#> 4: R1 D50 0.5 PSA 2779.9961
#> 5: R1 D50 1.0 PSA 780.4739
#> 6: R1 D50 2.0 PSA 129.8907widenPSL() projects the spectra to columns named
<ID>.<OCID> without changing values or
units.
widenPSL(PSL, by = "RecordID")[, .(Tn, PSA.H1, PSA.H2, PSA.D50)]
#> Tn PSA.H1 PSA.H2 PSA.D50
#> <num> <num> <num> <num>
#> 1: 0.0 1198.6235 903.18818 1344.4011
#> 2: 0.1 1231.6994 953.77134 1389.0485
#> 3: 0.2 1314.2224 1189.62599 1603.9814
#> 4: 0.5 2712.2604 1563.53961 2779.9961
#> 5: 1.0 956.5630 170.66027 780.4739
#> 6: 2.0 150.1647 37.27773 129.8907measureIntensitiesTSL() computes scalar intensity
measures for each record, component and signal domain after conversion
to the requested length basis; the units column identifies
the units of each measure.
IML <- measureIntensitiesTSL(Histories, by = "RecordID",
unitsSource = "mm", unitsTarget = "mm")
IML[IM %in% c("PGA", "AI", "D0595", "CAV", "PGV", "PGD")]
#> RecordID OCID ID IM value units
#> <char> <char> <char> <char> <num> <char>
#> 1: R1 H1 AT AI 867.247133 mm /s
#> 2: R1 H2 AT AI 490.393911 mm /s
#> 3: R1 H1 AT PGA 1199.068679 mm /s2
#> 4: R1 H2 AT PGA 899.163994 mm /s2
#> 5: R1 H1 AT D0595 9.400000 s
#> 6: R1 H2 AT D0595 9.375000 s
#> 7: R1 H1 AT CAV 7649.761708 mm /s
#> 8: R1 H2 AT CAV 5762.852857 mm /s
#> 9: R1 H1 VT PGV 130.016874 mm /s
#> 10: R1 H2 VT PGV 57.569217 mm /s
#> 11: R1 H1 DT PGD 14.212090 mm
#> 12: R1 H2 DT PGD 3.686879 mmProject functions read their JSON contract relative to
root and publish the products described in their help.
selectProject() filters index metadata and writes the
selected records and the candidate pool. Its contract requires
owners, path.selection and the index paths;
processID defaults to raw and
oneRecordPerEvent to TRUE, and PGA criteria require
units.PGA. The project below is written under
tempdir(): a processed index of four records from two
owners, and a contract with one criterion on magnitude and PGA.
Root <- file.path(tempdir(), "gmsp-introduction")
dir.create(file.path(Root, "database", "index"), recursive = TRUE)
dir.create(file.path(Root, "database", "process", "demo"), recursive = TRUE)
Index <- data.table(
RecordID = c("101", "102", "103", "104"), OwnerID = c("A", "B", "A", "B"),
EventID = c("2019a", "2019b", "2021a", "2022a"),
EventID.USGS = c("us01", "us01", "us02", "us03"),
StationID = c("S1", "S2", "S3", "S4"), DIR = "H1", OCID = "H1",
PGA = c(2000, 3500, 1500, 900), ProcessID = "demo",
EventMagnitude = c(6.8, 6.9, 7.1, 5.6),
EventMagnitudeType = c("Mw", "mww", "ML", "Mw")
)
fwrite(Index, file.path(Root, "database", "process", "demo", "MasterIndex.csv"))
jsonlite::write_json(list(
processID = "demo", owners = c("A", "B"),
units = list(PGA = "g"),
criteria = list(list(criterion = "strong",
EventMagnitude = c(6, 8), PGA = c(0.1, 0.5))),
path = list(index = "database/index", selection = "selection")
), file.path(Root, "select.json"), auto_unbox = TRUE, pretty = TRUE)The index stores PGA in mm/s^2, and units.PGA declares
the criterion bounds in g.
selectProject(file = "select.json", root = Root)
#> Selection published.
Selection <- fread(file.path(Root, "selection", "selection.csv"))
Candidates <- fread(file.path(Root, "selection", "candidates.csv"))
Selection[, .(RecordID, EventID.USGS, EventMagnitude, EventMagnitudeType, PGA)]
#> RecordID EventID.USGS EventMagnitude EventMagnitudeType PGA
#> <int> <char> <num> <char> <int>
#> 1: 102 us01 6.9 mww 3500
#> 2: 103 us02 7.1 Mw 1500
Candidates[, .(RecordID, EventID.USGS, PGA)]
#> RecordID EventID.USGS PGA
#> <int> <char> <int>
#> 1: 101 us01 2000
#> 2: 102 us01 3500
#> 3: 103 us02 1500Records belong to the same earthquake when their
EventID.USGS is equal, and with
oneRecordPerEvent TRUE selection.csv keeps one
record per earthquake, the one with the largest PGA, while
candidates.csv keeps every record that meets a criterion.
Before the criteria are applied, EventMagnitude becomes a
moment magnitude where one is published or an Ms is converted, and an ML
is taken as Mw without conversion. summary.json counts the
selected records and events with the same identity:
rankProject() scores candidate spectra against a target
band and writes a selection, processProject() publishes
record products and process tables from the stages gmsp,
trim, imf, imfTS and
psa, matchProject() composes modal shaping,
record scaling and population summaries, and
exportProject() exports project products. The articles Provider
files, raw records and indexes, Histories,
intensities and spectra, Selection,
spectral ranking and matching, Modal
decomposition and reconstruction and Project
steps, contracts and products describe the families of functions,
from the raw database to the run of a project. help(package = "gmsp")
lists the functions, and their help documents arguments, returned
objects and effects. The function reference and the books Mathematical
Basis and Numerical
Implementation are published at https://kanameishi.github.io/gmsp/.