support for planar loader

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. */