diff --git a/src/main/java/scratch/kevin/nshm23/prvi/SubductionDefModConvert.java b/src/main/java/scratch/kevin/nshm23/prvi/SubductionDefModConvert.java index 0917f5bb..76315f99 100644 --- a/src/main/java/scratch/kevin/nshm23/prvi/SubductionDefModConvert.java +++ b/src/main/java/scratch/kevin/nshm23/prvi/SubductionDefModConvert.java @@ -13,12 +13,16 @@ import java.util.Comparator; import java.util.HashMap; import java.util.HashSet; +import java.util.LinkedHashMap; import java.util.List; import java.util.Map; import org.jfree.chart.annotations.XYTextAnnotation; import org.jfree.chart.ui.TextAnchor; +import org.jfree.data.Range; import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc; +import org.opensha.commons.data.function.DefaultXY_DataSet; +import org.opensha.commons.data.function.XY_DataSet; import org.opensha.commons.geo.Location; import org.opensha.commons.geo.LocationList; import org.opensha.commons.geo.LocationUtils; @@ -29,13 +33,20 @@ import org.opensha.commons.geo.json.Geometry.DepthSerializationType; import org.opensha.commons.geo.json.Geometry.MultiLineString; import org.opensha.commons.gui.plot.GeographicMapMaker; +import org.opensha.commons.gui.plot.HeadlessGraphPanel; +import org.opensha.commons.gui.plot.PlotCurveCharacterstics; +import org.opensha.commons.gui.plot.PlotLineType; +import org.opensha.commons.gui.plot.PlotSpec; import org.opensha.commons.gui.plot.PlotSymbol; +import org.opensha.commons.gui.plot.PlotUtils; import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; import org.opensha.commons.util.FaultUtils; import org.opensha.commons.util.MarkdownUtils; import org.opensha.commons.util.MarkdownUtils.TableBuilder; import org.opensha.commons.util.cpt.CPT; import org.opensha.commons.util.DataUtils.MinMaxAveTracker; +import org.opensha.sha.earthquake.faultSysSolution.util.SubSectionBuilder; +import org.opensha.sha.earthquake.faultSysSolution.util.minisections.MinisectionMappings; import org.opensha.sha.earthquake.faultSysSolution.util.minisections.MinisectionSlipRecord; import org.opensha.sha.faultSurface.FaultSection; import org.opensha.sha.faultSurface.FaultTrace; @@ -44,6 +55,7 @@ import com.google.common.base.Preconditions; import com.google.common.collect.HashBasedTable; import com.google.common.collect.Table; +import com.google.common.primitives.Ints; import com.google.gson.Gson; import com.google.gson.JsonIOException; import com.google.gson.JsonSyntaxException; @@ -64,6 +76,7 @@ public static void main(String[] args) throws IOException { boolean[] largePolyBools = { true, false }; boolean interpolate = true; + boolean interpSymmetry = true; File debugDir = new File("/tmp/sub_fm_dm_debug"); Preconditions.checkState(debugDir.exists() || debugDir.mkdir()); @@ -123,6 +136,14 @@ public static void main(String[] args) throws IOException { stitchInIDs.put(name2, inIDs2); stitchOutIDs.put(name2, outID2); + Map interpInIDs = null; + + if (interpolate) { + interpInIDs = new HashMap<>(); + interpInIDs.put("Hispaniola - PRVI - Lesser Antilles", Ints.concat(inIDs0, inIDs1)); + interpInIDs.put("Muertos", inIDs2); + } + Map stitchFeatureSets = new HashMap<>(); for (String name : stitchInIDs.keySet()) stitchFeatureSets.put(name, new Feature[stitchInIDs.get(name).length]); @@ -147,6 +168,7 @@ public static void main(String[] args) throws IOException { List outFeatures = new ArrayList<>(); Map> minisectionRecsMap = new HashMap<>(); + Map> minisectionRecsMapToOrigIDs = new HashMap<>(); for (String name : stitchInIDs.keySet()) { int outID = stitchOutIDs.get(name); @@ -161,11 +183,11 @@ public static void main(String[] args) throws IOException { int prevID = -1; - List minisectionDASs = new ArrayList<>(); - ArbitrarilyDiscretizedFunc dasSlipFunc = new ArbitrarilyDiscretizedFunc(); - ArbitrarilyDiscretizedFunc dasSlipUncertFunc = new ArbitrarilyDiscretizedFunc(); - ArbitrarilyDiscretizedFunc dasRakeFunc = new ArbitrarilyDiscretizedFunc(); - double runningDAS = 0d; +// List minisectionDASs = new ArrayList<>(); +// ArbitrarilyDiscretizedFunc dasSlipFunc = new ArbitrarilyDiscretizedFunc(); +// ArbitrarilyDiscretizedFunc dasSlipUncertFunc = new ArbitrarilyDiscretizedFunc(); +// ArbitrarilyDiscretizedFunc dasRakeFunc = new ArbitrarilyDiscretizedFunc(); +// double runningDAS = 0d; for (Feature feature : features) { Preconditions.checkState(feature.geometry.type == GeoJSON_Type.MultiLineString); @@ -235,43 +257,22 @@ public static void main(String[] args) throws IOException { double slipUncert = props.getDouble(rateUncertPropName, Double.NaN); double rake = props.getDouble(rakePropName, Double.NaN); - if (runningDAS == 0d) { - dasSlipFunc.set(0d, slip); - dasSlipUncertFunc.set(0d, slipUncert); - dasRakeFunc.set(0d, rake); - } - - double myStartDAS = runningDAS; + List myMinis = new ArrayList<>(upper.size()-1); for (int i=1; i= 0 && das <= runningDAS); - minisectionRecs.set(i, new MinisectionSlipRecord(rec.parentID, rec.minisectionID, - rec.startLoc, rec.endLoc, - dasRakeFunc.getInterpolatedY(das), - dasSlipFunc.getInterpolatedY(das), - dasSlipUncertFunc.getInterpolatedY(das))); + minisectionRecsMap.put(outID, minisectionRecs); + } + + if (interpolate) { + System.out.println("Interpolating"); + + // build original subsections and calculate original moment rates + List fullSects = new ArrayList<>(); + for (Feature feature : outFeatures) + fullSects.add(GeoJSONFaultSection.fromFeature(feature)); + List subSects = SubSectionBuilder.buildSubSects(fullSects, 2, 0.5, 30); + MinisectionMappings mappings = new MinisectionMappings(fullSects, subSects); + mappings.mapDefModelMinisToSubSects(minisectionRecsMap); + Map origMoRates = new HashMap<>(); + for (FaultSection sect : subSects) { + String name = sect.getParentSectionName(); + Double moRate = origMoRates.get(name); + if (moRate == null) + moRate = 0d; + origMoRates.put(name, moRate+sect.calcMomentRate(false)); + } + + for (String name : interpInIDs.keySet()) { + int[] ids = interpInIDs.get(name); + System.out.println("Interpolating "+name); + + // figure out section bounds in terms of DAS + double runningDAS = 0d; + List interpMinis = new ArrayList<>(); + List miniMidDASs = new ArrayList<>(); + + List sectStartDASs = new ArrayList<>(); + List sectMidDASs = new ArrayList<>(); + List sectEndDASs = new ArrayList<>(); + List sectRefMinis = new ArrayList<>(); + + for (int id : ids) { + List minis = minisectionRecsMapToOrigIDs.get(id); + Preconditions.checkNotNull(minis); + + double start = runningDAS; + + for (MinisectionSlipRecord mini : minis) { + double len = LocationUtils.linearDistanceFast(mini.startLoc, mini.endLoc); + runningDAS += len; + + interpMinis.add(mini); + miniMidDASs.add(runningDAS - 0.5*len); + } + double end = runningDAS; + double middle = 0.5*(start + end); + runningDAS = end; + + sectRefMinis.add(minis.get(0)); + sectStartDASs.add(start); + sectMidDASs.add(middle); + sectEndDASs.add(end); + } + + System.out.println("Total DAS for "+name+": "+runningDAS); + + // build slip/rake as a function of DAS + ArbitrarilyDiscretizedFunc dasSlipFunc = new ArbitrarilyDiscretizedFunc(); + ArbitrarilyDiscretizedFunc dasSlipUncertFunc = new ArbitrarilyDiscretizedFunc(); + ArbitrarilyDiscretizedFunc dasRakeFunc = new ArbitrarilyDiscretizedFunc(); + + for (int i=0; i prevHalfLen) { + // my portion is longer, truncate + double interpX = start + prevHalfLen; + dasSlipFunc.set(interpX, mini.slipRate); + dasSlipUncertFunc.set(interpX, mini.slipRateStdDev); + dasRakeFunc.set(interpX, mini.rake); + } else { + // previous portion was longer, truncate + MinisectionSlipRecord prev = sectRefMinis.get(i-1); + double interpX = start - myHalfLen; + dasSlipFunc.set(interpX, prev.slipRate); + dasSlipUncertFunc.set(interpX, prev.slipRateStdDev); + dasRakeFunc.set(interpX, prev.rake); + } + } + + // add the middle + dasSlipFunc.set(middle, mini.slipRate); + dasSlipUncertFunc.set(middle, mini.slipRateStdDev); + dasRakeFunc.set(middle, mini.rake); + + if (i == sectRefMinis.size()-1) { + // last, extend to the end + dasSlipFunc.set(end, mini.slipRate); + dasSlipUncertFunc.set(end, mini.slipRateStdDev); + dasRakeFunc.set(end, mini.rake); + } + } + + // apply interpolation + for (int i=0; i origRecs = minisectionRecsMap.get(mini.parentID); + Preconditions.checkState(mini.minisectionID < origRecs.size()); + MinisectionSlipRecord orig = origRecs.get(mini.minisectionID); + Preconditions.checkState(orig.equals(mini)); + origRecs.set(mini.minisectionID, modMini); } + + // plot them + String interpPrefix = prefix+"_interp_"+name.replaceAll("\\W+", "_"); + debugWriteInterpFuncsOrders(debugDir, interpPrefix, name, + dasSlipFunc, dasSlipUncertFunc, dasRakeFunc, sectStartDASs); } - minisectionRecsMap.put(outID, minisectionRecs); + // calculate modified moment rates + mappings.mapDefModelMinisToSubSects(minisectionRecsMap); + Map modMoRates = new LinkedHashMap<>(); + for (FaultSection sect : subSects) { + String name = sect.getParentSectionName(); + Double moRate = modMoRates.get(name); + if (moRate == null) + moRate = 0d; + modMoRates.put(name, moRate+sect.calcMomentRate(false)); + } + + System.out.println("Interpolation Mo-Rate changes for "+title); + for (String name : modMoRates.keySet()) { + double orig = origMoRates.get(name); + double mod = modMoRates.get(name); + System.out.println("\t"+name+":\t"+moDF.format(orig)+" -> "+moDF.format(mod)+" ("+pDF.format((mod-orig)/orig)+")"); + } } if (fmOutputFile != null) @@ -360,6 +500,9 @@ public static void main(String[] args) throws IOException { MarkdownUtils.writeReadmeAndHTML(lines, debugDir); } + private static final DecimalFormat moDF = new DecimalFormat("0.00E0"); + private static final DecimalFormat pDF = new DecimalFormat("0.00%"); + private static String depthStr(FaultTrace trace) { MinMaxAveTracker track = new MinMaxAveTracker(); for (Location loc : trace) @@ -450,5 +593,62 @@ private static void debugWriteSectOrders(List features, File outputDir, mapMaker.plotSectScalars(rakes, rakeCPT, "Rake Angle"); mapMaker.plot(outputDir, prefix+"_rakes", title); } + + private static void debugWriteInterpFuncsOrders(File outputDir, String prefix, String name, + ArbitrarilyDiscretizedFunc dasSlipFunc, ArbitrarilyDiscretizedFunc dasSlipUncertFunc, + ArbitrarilyDiscretizedFunc dasRakeFunc, List sectStartDASs) throws IOException { + + List specs = new ArrayList<>(3); + List yRanges = new ArrayList<>(3); + + PlotCurveCharacterstics startChar = new PlotCurveCharacterstics(PlotLineType.DASHED, 1f, Color.GRAY); + + for (int i=0; i<3; i++) { + ArbitrarilyDiscretizedFunc inFunc; + String yAxisLabel; + switch (i) { + case 0: + inFunc = dasSlipFunc; + yAxisLabel = "Slip Rate (mm/yr)"; + break; + case 1: + inFunc = dasSlipUncertFunc; + yAxisLabel = "Slip Rate Uncertainty (mm/yr)"; + break; + case 2: + inFunc = dasRakeFunc; + yAxisLabel = "Rake Angle"; + break; + default: + throw new IllegalStateException(); + } + + double maxY = inFunc.getMaxY()*1.05; + yRanges.add(new Range(0d, maxY)); + + List funcs = new ArrayList<>(); + List chars = new ArrayList<>(); + + for (int s=1; s