From 5e383283411e8103cf73950a2b86b966c41f8adc Mon Sep 17 00:00:00 2001
From: Peter Powers <pmpowers@usgs.gov>
Date: Fri, 28 Mar 2025 14:59:11 -0600
Subject: [PATCH] support for planar loader
---
.../earthquake/nshmp/model/GridLoader.java | 198 +++++++++++++++++-
1 file changed, 194 insertions(+), 4 deletions(-)
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java b/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java
index 11a1a09d..ca45dd3f 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java
@@ -4,6 +4,7 @@ import static com.google.common.base.Preconditions.checkState;
import static gov.usgs.earthquake.nshmp.model.Deserialize.MFD_TREE;
import static java.lang.Math.log10;
import static java.lang.Math.pow;
+import static java.util.stream.Collectors.toList;
import static java.util.stream.Collectors.toUnmodifiableList;
import java.nio.file.Path;
@@ -65,8 +66,13 @@ class GridLoader {
static final String LON = "lon";
static final String LAT = "lat";
static final String DEPTH = "depth";
- static final String STRIKE = "strike";
+ static final String MAG = "m";
static final String RATE = "rate";
+ static final String STRIKE = "strike";
+ static final String DIP = "dip";
+ static final String UPPER_DEPTH = "zTor";
+ static final String LOWER_DEPTH = "zBor";
+ static final String LENGTH = "length";
static final String NAME = "name";
static final String PDF = "pdf";
static final String MMAX = "mMax";
@@ -82,8 +88,13 @@ class GridLoader {
String mfdTreeKey,
ModelData data) {
- SlabLoader loader = new SlabLoader(id, name, featureId, mfdTreeKey, data);
- return loader.createSlabRuptureSet();
+ if (data.spatialPdf()) {
+ SlabLoader loader = new SlabLoader(id, name, featureId, mfdTreeKey, data);
+ return loader.createSlabRuptureSet();
+ }
+ PlanarLoader loader = new PlanarLoader(id, name, featureId, mfdTreeKey, data);
+ return loader.createRuptureSet();
+
}
static GridRuptureSet createGridRuptureSet(
@@ -107,6 +118,183 @@ class GridLoader {
return loader.createRuptureSet();
}
+ /*
+ * Explicit INCR MFD with one line per magnitude and necessary data to realize
+ * each rupture as a finite, planar surface. Introduced for PRVI NSHM for
+ * subduction interface gridded seismicity.
+ */
+ private static class PlanarLoader {
+
+ private final int id;
+ private final String name;
+ private final ModelData data;
+ private final SourceFeature.Grid feature;
+ private final Path mfdsPath;
+ private final String mfdTreeKey;
+ private final FeatureData featureData;
+
+ private PlanarLoader(
+ int id,
+ String name,
+ int featureId,
+ String mfdTreeKey,
+ ModelData data) {
+
+ this.id = id;
+ this.name = name;
+ this.data = data;
+ this.feature = data.gridFeatureMap().orElseThrow().get(featureId);
+ this.mfdsPath = data.gridRateFile().orElseThrow();
+ this.mfdTreeKey = mfdTreeKey;
+ this.featureData = new FeatureData();
+ processNodes();
+ }
+
+ private SlabRuptureSet createRuptureSet() {
+ return SlabRuptureSet.builder()
+ .name(name)
+ .id(id)
+ .modelData(data)
+ .feature(feature)
+ .nodeData(
+ featureData.locations,
+ featureData.mfds,
+ featureData.mfdsTree,
+ featureData.ruptureDataLists)
+ .build();
+ }
+
+ private void processNodes() {
+
+ DelimitedData ratefile = DelimitedData.comma(mfdsPath);
+
+ /* Because we have explicit MFDs, we expect a single branch mfd-tree. */
+ LogicTree<Mfd.Properties> mfdTree = data.mfdMap().orElseThrow().get(mfdTreeKey);
+ checkState(mfdTree.size() == 1);
+ Mfd.Properties mfdProps = mfdTree.get(0).value();
+
+ List<Record> records;
+ try (Stream<Record> s = ratefile.records()) {
+ records = s.collect(toList());
+ }
+ List<List<Record>> nodeRecordsList = splitRecords(records);
+ // process record group for each grid node
+ for (List<Record> nodeRecords : nodeRecordsList) {
+ processRecords(nodeRecords, mfdProps);
+ }
+
+ featureData.mfdsTree = LogicTree.singleton(
+ feature.name,
+ mfdTree.get(0).id(),
+ featureData.mfds);
+ }
+
+ /*
+ * Break records into groups by node; method assumes no repeating lat-lons
+ * with different depths
+ */
+ private List<List<Record>> splitRecords(List<Record> records) {
+ List<List<Record>> nodeRecordLists = new ArrayList<>();
+ List<Record> nodeRecords = new ArrayList<>();
+ for (int i = 0; i < records.size(); i++) {
+ Record cRecord = records.get(i);
+ String cLon = cRecord.get(Key.LON);
+ String cLat = cRecord.get(Key.LAT);
+ if (i == records.size() - 1) {
+ nodeRecords.add(cRecord);
+ nodeRecordLists.add(nodeRecords);
+ } else {
+ nodeRecords.add(cRecord);
+ // check next record and possibly reset list
+ Record nRecord = records.get(i + 1);
+ String nLon = nRecord.get(Key.LON);
+ String nLat = nRecord.get(Key.LAT);
+ if (!nLon.equals(cLon) || !nLat.equals(cLat)) {
+ nodeRecordLists.add(nodeRecords);
+ nodeRecords = new ArrayList<>();
+ }
+ }
+ }
+ return nodeRecordLists;
+ }
+
+ private void processRecords(List<Record> records, Mfd.Properties incrMfdProps) {
+
+ Location loc = Location.create(
+ records.get(0).getDouble(Key.LON),
+ records.get(0).getDouble(Key.LAT),
+ records.get(0).getDouble(Key.DEPTH));
+
+ if (!feature.contains(loc)) {
+ // TODO clean
+ System.out.println("Point outside poly");
+ return;
+ }
+ featureData.locations.add(loc);
+
+ Mfd.Builder mfdBuilder = incrMfdProps.toBuilder();
+ IncrRates incrRates = new IncrRates(records);
+ mfdBuilder.transform(incrRates::set);
+ featureData.mfds.add(mfdBuilder.build());
+
+ List<RuptureData> rupDataList = new ArrayList<>();
+ for (int i = 0; i < records.size(); i++) {
+ Record r = records.get(i);
+ RuptureData rupDat = new RuptureData(
+ r.getDouble(Key.STRIKE),
+ r.getDouble(Key.DIP),
+ r.getDouble(Key.UPPER_DEPTH),
+ r.getDouble(Key.LOWER_DEPTH),
+ r.getDouble(Key.LENGTH));
+ rupDataList.add(rupDat);
+ }
+ featureData.ruptureDataLists.add(rupDataList);
+ }
+
+ private static class IncrRates {
+
+ final List<Record> records;
+
+ IncrRates(List<Record> records) {
+ this.records = records;
+ }
+
+ /*
+ * Using XyPoint.index() makes the assumption that there is agreement in
+ * magnitude range and ascending order. The number of records for each
+ * node can be less than the
+ */
+ //
+ void set(XyPoint mBin) {
+ int index = mBin.index();
+ if (index >= records.size()) {
+ return;
+ }
+ // Record record = records.get(index);
+ mBin.set(records.get(index).getDouble(Key.RATE));
+ double m = records.get(index).getDouble(Key.MAG);
+ checkState(m == mBin.x(), "Magnitudes don't match; mfd: %s record: %s", mBin.x(), m);
+ }
+ }
+ }
+
+ static class RuptureData {
+
+ final double strike;
+ final double dip;
+ final double zTor;
+ final double zBor;
+ final double length;
+
+ RuptureData(double strike, double dip, double zTor, double zBor, double length) {
+ this.strike = strike;
+ this.dip = dip;
+ this.zTor = zTor;
+ this.zBor = zBor;
+ this.length = length;
+ }
+ }
+
/* Explicit INCR MFD (fault systems, e.g. UCERF3). */
private static class IncrLoader {
@@ -350,7 +538,8 @@ class GridLoader {
.nodeData(
featureData.locations,
featureData.mfds,
- featureData.mfdsTree);
+ featureData.mfdsTree,
+ featureData.ruptureDataLists);
return b.build();
}
@@ -443,6 +632,7 @@ class GridLoader {
/* transposed and combined data */
LogicTree<List<Mfd>> mfdsTree;
List<Mfd> mfds = new ArrayList<>();
+ List<List<RuptureData>> ruptureDataLists = new ArrayList<>();
}
/* Truncate special GR cases; e.g. WUS double counting. */
--
GitLab