diff --git a/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java b/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java index c84c187f..4e3c35cd 100644 --- a/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java +++ b/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java @@ -625,35 +625,35 @@ public static void main(String[] args) throws IOException { * PRVI25 logic tree * TODO (this is a just a marker to find this part quickly, not an actual todo) */ -// List> levels = PRVI25_LogicTreeBranch.levelsOnFault; -// dirName += "-prvi25_crustal_branches"; -// double avgNumRups = 50000; -// gmpes = new AttenRelRef[] { AttenRelRef.USGS_PRVI_ACTIVE }; -// -// // random DM sampling -// levels = new ArrayList<>(levels); -// int origNumLevels = levels.size(); -// for (int i=levels.size(); --i>=0;) -// if (levels.get(i).getNodes().get(0) instanceof PRVI25_CrustalDeformationModels) -// levels.remove(i); -// Preconditions.checkState(levels.size() == origNumLevels -1); -// individualRandomLevels.add(new PRVI25_CrustalRandomlySampledDeformationModelLevel()); -// samplingBranchCountMultiplier = 5; // 5 for each branch -// dirName += "-dmSample"; -// if (samplingBranchCountMultiplier > 1) -// dirName += samplingBranchCountMultiplier+"x"; + List> levels = PRVI25_LogicTreeBranch.levelsOnFault; + dirName += "-prvi25_crustal_branches"; + double avgNumRups = 50000; + gmpes = new AttenRelRef[] { AttenRelRef.USGS_PRVI_ACTIVE }; + + // random DM sampling + levels = new ArrayList<>(levels); + int origNumLevels = levels.size(); + for (int i=levels.size(); --i>=0;) + if (levels.get(i).getNodes().get(0) instanceof PRVI25_CrustalDeformationModels) + levels.remove(i); + Preconditions.checkState(levels.size() == origNumLevels -1); + individualRandomLevels.add(new PRVI25_CrustalRandomlySampledDeformationModelLevel()); + samplingBranchCountMultiplier = 5; // 5 for each branch + dirName += "-dmSample"; + if (samplingBranchCountMultiplier > 1) + dirName += samplingBranchCountMultiplier+"x"; - List> levels = PRVI25_LogicTreeBranch.levelsSubduction; - dirName += "-prvi25_subduction_branches"; - double avgNumRups = 10000; - gmpes = new AttenRelRef[] { AttenRelRef.USGS_PRVI_INTERFACE, AttenRelRef.USGS_PRVI_SLAB }; +// List> levels = PRVI25_LogicTreeBranch.levelsSubduction; +// dirName += "-prvi25_subduction_branches"; +// double avgNumRups = 10000; +// gmpes = new AttenRelRef[] { AttenRelRef.USGS_PRVI_INTERFACE, AttenRelRef.USGS_PRVI_SLAB }; // levels = new ArrayList<>(levels); // levels.add(NSHM23_LogicTreeBranch.SUB_SECT_CONSTR); // dirName += "-proxyGriddedTests"; -// Class factoryClass = PRVI25_InvConfigFactory.class; + Class factoryClass = PRVI25_InvConfigFactory.class; // Class factoryClass = PRVI25_InvConfigFactory.GriddedUseM1Bounds.class; // dirName += "-grid_bounds_m1"; @@ -670,8 +670,8 @@ public static void main(String[] args) throws IOException { // Class factoryClass = PRVI25_InvConfigFactory.RateBalanceAndLimitCrustalBelowObserved_0p9.class; // dirName += "-limit_below_obs_constraint-grided_rate_balancing"; - Class factoryClass = PRVI25_InvConfigFactory.GriddedForceSlab2Depths.class; - dirName += "-gridded_use_slab2"; +// Class factoryClass = PRVI25_InvConfigFactory.GriddedForceSlab2Depths.class; +// dirName += "-gridded_use_slab2"; if (!factoryClass.equals(PRVI25_InvConfigFactory.class)) { // try instantiate it to make sure we get any static modifiers that might change branch weights diff --git a/src/main/java/scratch/kevin/prvi25/FaultSystemLineIntegralCalculator.java b/src/main/java/scratch/kevin/prvi25/FaultSystemLineIntegralCalculator.java index 121b40aa..43fe6774 100644 --- a/src/main/java/scratch/kevin/prvi25/FaultSystemLineIntegralCalculator.java +++ b/src/main/java/scratch/kevin/prvi25/FaultSystemLineIntegralCalculator.java @@ -483,10 +483,12 @@ public GeographicMapMaker buildMapPlot(boolean labelMagnitudes, boolean plotVect if (result == integrals.get(0)) { double integralMiddleLat = 0.5*(result.startLoc.lat + result.endLoc.lat); double integralMiddleLon = 0.5*(result.startLoc.lon + result.endLoc.lon); + double furthestLeftLon = Math.max(lowerLeft.lon+0.25, lowerLeft.lon + 0.05*(upperRight.lon - lowerLeft.lon)); + System.out.println("Furthest left="+furthestLeftLon+", integral left="+integralMiddleLon); XYTextAnnotation ann = new XYTextAnnotation("Summation Direction", integralMiddleLon, integralMiddleLat); ann.setFont(textAnnFont); - ann.setTextAnchor(TextAnchor.TOP_CENTER); - ann.setRotationAnchor(TextAnchor.TOP_CENTER); + ann.setTextAnchor(TextAnchor.BOTTOM_CENTER); + ann.setRotationAnchor(TextAnchor.BOTTOM_CENTER); double vectorAngle = result.azimuth; while (vectorAngle > 180d) vectorAngle -= 360d; @@ -496,6 +498,13 @@ public GeographicMapMaker buildMapPlot(boolean labelMagnitudes, boolean plotVect ann.setTextAnchor(TextAnchor.BOTTOM_CENTER); ann.setRotationAnchor(TextAnchor.BOTTOM_CENTER); } +// } else if (integralMiddleLon < furthestLeftLon) { +//// } else { +// // too far to the left, put it on the other side +//// vectorAngle -= 180; +// ann.setTextAnchor(TextAnchor.BOTTOM_CENTER); +// ann.setRotationAnchor(TextAnchor.BOTTOM_CENTER); +// } double angleRad = Math.toRadians(vectorAngle + 90); ann.setRotationAngle(mapMaker.getRotationAngleCorrectedForAspectRatio(angleRad)); mapMaker.addAnnotation(ann); @@ -587,7 +596,7 @@ public GeographicMapMaker buildMapPlot(boolean labelMagnitudes, boolean plotVect arrowLen = 0.2*maxVectDist; arrowThickness = 4f; } else { - arrowLen = 0.2*maxVectDist; + arrowLen = 0.15*maxVectDist; arrowThickness = 2f; } mapMaker.setFillArrowheads(true); diff --git a/src/main/java/scratch/kevin/prvi25/PRVI_SubductionTestRupSetBuilder.java b/src/main/java/scratch/kevin/prvi25/PRVI_SubductionTestRupSetBuilder.java index 3901a6e5..85630ae7 100644 --- a/src/main/java/scratch/kevin/prvi25/PRVI_SubductionTestRupSetBuilder.java +++ b/src/main/java/scratch/kevin/prvi25/PRVI_SubductionTestRupSetBuilder.java @@ -33,7 +33,8 @@ public static void main(String[] args) throws IOException { // PRVI25_SubductionFaultModels fm = PRVI25_SubductionFaultModels.PRVI_SUB_FM_SMALL; PRVI25_SubductionDeformationModels dm = PRVI25_SubductionDeformationModels.FULL; - File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_Subduction/figures/fault_model"); +// File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_Subduction/figures/fault_model"); + File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_ERF/prvi25-erf-paper/Figures/sub_fm"); Preconditions.checkState(outputDir.exists() || outputDir.mkdir()); // first just write out the fault model diff --git a/src/main/java/scratch/kevin/prvi25/SubductionDefModConvert.java b/src/main/java/scratch/kevin/prvi25/SubductionDefModConvert.java index 7811f1c3..287d2877 100644 --- a/src/main/java/scratch/kevin/prvi25/SubductionDefModConvert.java +++ b/src/main/java/scratch/kevin/prvi25/SubductionDefModConvert.java @@ -95,8 +95,9 @@ public static void main(String[] args) throws IOException { boolean interpolate = true; boolean interpSymmetry = true; - File debugDir = new File("/home/kevin/Documents/papers/2024_PRVI_Subduction/figures/def_model"); - Preconditions.checkState(debugDir.exists() || debugDir.mkdir()); +// File plotOutputDir = new File("/home/kevin/Documents/papers/2024_PRVI_Subduction/figures/def_model"); + File plotOutputDir = new File("/home/kevin/Documents/papers/2024_PRVI_ERF/prvi25-erf-paper/Figures/sub_dm"); + Preconditions.checkState(plotOutputDir.exists() || plotOutputDir.mkdir()); Table debugPrefixes = HashBasedTable.create(); @@ -138,7 +139,7 @@ public static void main(String[] args) throws IOException { String fmTitle = geomTitle; String fmDmTitle = geomTitle+", "+rateTitle; debugPrefixes.put(fullRate, largePolys, prefix); - debugWriteSectOrders(inFeatures, debugDir, prefix, fmTitle, fmDmTitle, + debugWriteSectOrders(inFeatures, plotOutputDir, prefix, fmTitle, fmDmTitle, ratePropName, rateUncertPropName, rakePropName); Map stitchInIDs = new HashMap<>(); @@ -447,7 +448,7 @@ public static void main(String[] args) throws IOException { // plot them String interpPrefix = prefix+"_interp_"+name.replaceAll("\\W+", "_"); - debugWriteInterpFuncsOrders(debugDir, interpPrefix, name, + debugWriteInterpFuncsOrders(plotOutputDir, interpPrefix, name, dasSlipFunc, dasSlipUncertFunc, dasRakeFunc, sectStartDASs, sectMidDASs, interpRanges); } @@ -531,7 +532,7 @@ public static void main(String[] args) throws IOException { lines.add(tocIndex, "## Table Of Contents"); // write markdown - MarkdownUtils.writeReadmeAndHTML(lines, debugDir); + MarkdownUtils.writeReadmeAndHTML(lines, plotOutputDir); } private static final DecimalFormat moDF = new DecimalFormat("0.00E0"); @@ -652,8 +653,13 @@ private static void debugWriteSectOrders(List features, File outputDir, for (int s=0; s clustersUnsorted = rupSet.requireModule(ConnectivityClusters.class).get(); + + // see if any clusters rupture completely + System.out.println("Looking for clusters that rupture fully..."); + boolean anyFound = false; + for (ConnectivityCluster cluster : clustersUnsorted) { + HashSet parents = new HashSet<>(cluster.getParentSectIDs()); + if (parents.size() == 1) + continue; + Set sects = cluster.getSectIDs(); + for (List rupSects : rupSet.getSectionIndicesForAllRups()) { + if (rupSects.size() == sects.size() && rupSects.containsAll(sects)) { + System.out.println("Cluster ruptures fully: "+sects); + anyFound = true; + System.out.println("Parents:"); + for (int sectID: rupSects) { + FaultSection sect = rupSet.getFaultSectionData(sectID); + if (parents.contains(sect.getParentSectionId())) { + System.out.println("\t"+sect.getParentSectionId()+". "+sect.getParentSectionName()); + parents.remove(sect.getParentSectionId()); + } + } + } + } + } + if (!anyFound) + System.out.println("No clusters rupture fully"); + + GeographicMapMaker plotter = FaultSectionConnectionsPlot.buildConnectedClustersPlot(rupSet, sol, + SlipRateFigures.CRUSTAL_FAULT_MAP_REG, null, clustersUnsorted); + + plotter.setWritePDFs(true); + plotter.setWriteGeoJSON(false); + + plotter.plot(outputDir, "connectivity_clusters", " "); + + CrustalFaultNamesFigure.addStandardFaultLabels(plotter, rupSet.getFaultSectionDataList()); + + plotter.plot(outputDir, "connectivity_clusters_with_fault_names", " "); + } + +} diff --git a/src/main/java/scratch/kevin/prvi25/figures/CrustalFaultNamesFigure.java b/src/main/java/scratch/kevin/prvi25/figures/CrustalFaultNamesFigure.java new file mode 100644 index 00000000..b53da4af --- /dev/null +++ b/src/main/java/scratch/kevin/prvi25/figures/CrustalFaultNamesFigure.java @@ -0,0 +1,145 @@ +package scratch.kevin.prvi25.figures; + +import java.awt.Color; +import java.awt.Font; +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; +import java.util.Collection; +import java.util.HashMap; +import java.util.List; +import java.util.Map; + +import org.jfree.chart.annotations.XYTextAnnotation; +import org.jfree.chart.ui.TextAnchor; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.LocationList; +import org.opensha.commons.geo.LocationUtils; +import org.opensha.commons.gui.plot.GeographicMapMaker; +import org.opensha.commons.gui.plot.PlotCurveCharacterstics; +import org.opensha.commons.gui.plot.PlotLineType; +import org.opensha.commons.util.FaultUtils; +import org.opensha.sha.earthquake.rupForecastImpl.prvi25.logicTree.PRVI25_CrustalFaultModels; +import org.opensha.sha.faultSurface.FaultSection; +import org.opensha.sha.faultSurface.FaultTrace; + +import net.mahdilamb.colormap.Colors; + +class CrustalFaultNamesFigure { + + public static void main(String[] args) throws IOException { + File ouutputDir = new File("/home/kevin/Documents/papers/2024_PRVI_ERF/prvi25-erf-paper/Figures/crustal_fm"); + + PRVI25_CrustalFaultModels fm = PRVI25_CrustalFaultModels.PRVI_CRUSTAL_FM_V1p1; + +// List sects = fm.getFaultSections(); + List sects = fm.getDefaultDeformationModel().build(fm); + + GeographicMapMaker mapMaker = new GeographicMapMaker(SlipRateFigures.CRUSTAL_FAULT_MAP_REG, sects); + mapMaker.setWriteGeoJSON(false); + + addStandardFaultLabels(mapMaker, sects); + + mapMaker.setSectTraceChar(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, Colors.tab_blue)); + + mapMaker.plot(ouutputDir, "crustal_sect_names", " "); + } + + static void addStandardFaultLabels(GeographicMapMaker mapMaker, List sects) { + List anns = new ArrayList<>(); + List arrows = new ArrayList<>(); + Font labelFont = new Font(Font.SANS_SERIF, Font.BOLD, 14); +// buildLabelsAndArrows(sects, "Bunce", false, new Location(20, -66), +// labelFont, TextAnchor.BASELINE_CENTER, 0d, anns, arrows); + buildLabelsAndArrows(sects, List.of("Bunce 1", "Bunce 2", "Bunce 3", "Bunce 4", "Bunce 5"), "Bunce 1-5", false, new Location(19.75, -67.7), + labelFont, TextAnchor.BASELINE_CENTER, -8d, anns, arrows); + buildLabelsAndArrows(sects, List.of("Bunce 6", "Bunce 7", "Bunce 8"), "Bunce 6-8", false, new Location(19.45, -64), + labelFont, TextAnchor.TOP_CENTER, 8d, anns, arrows); + buildLabelsAndArrows(sects, "Anegada Passage", false, new Location(18.7, -64.7), + labelFont, TextAnchor.BASELINE_CENTER, -27, anns, arrows); + buildLabelsAndArrows(sects, "Main Ridge", false, new Location(19.2, -66.3), + labelFont, TextAnchor.TOP_CENTER, 0d, anns, arrows); + buildLabelsAndArrows(sects, "Septentrional", false, new Location(19.0, -69.5), + labelFont, TextAnchor.TOP_CENTER, 9d, anns, arrows); + buildLabelsAndArrows(sects, "South Lajas", false, new Location(17.5, -67.35), + labelFont, TextAnchor.TOP_CENTER, 0d, anns, arrows); + + mapMaker.addAnnotations(anns); + mapMaker.plotArrows(arrows, 6d, new Color(0, 0, 0, 180), 1f); + mapMaker.setFillArrowheads(true); + } + + static void buildLabelsAndArrows(List sects, String name, boolean stitch, + Location textLoc, Font font, TextAnchor textAnchor, double rotationAngle, + List anns, List arrows) { + buildLabelsAndArrows(sects, List.of(name), name, stitch, textLoc, font, textAnchor, rotationAngle, anns, arrows); + } + + static void buildLabelsAndArrows(List sects, List searchNames, String plotName, boolean stitch, + Location textLoc, Font font, TextAnchor textAnchor, double rotationAngle, + List anns, List arrows) { + Map> byParent = stitch ? null : new HashMap<>(); + List allCombined = stitch ? new ArrayList<>() : null; + for (FaultSection sect : sects) { + for (String name : searchNames) { + if (sect.getName().contains(name)) { + if (stitch) { + allCombined.add(sect); + } else { + int parent = sect.getParentSectionId(); + if (parent == -1) + // already parents + parent = sect.getSectionId(); + if (!byParent.containsKey(parent)) + byParent.put(parent, new ArrayList<>()); + byParent.get(parent).add(sect); + } + break; + } + } + } + + XYTextAnnotation ann = new XYTextAnnotation(plotName, textLoc.getLongitude(), textLoc.getLatitude()); + ann.setTextAnchor(textAnchor); + ann.setFont(font); + if (rotationAngle != 0d) + ann.setRotationAngle(GeographicMapMaker.getRotationAngleCorrectedForAspectRatio(Math.toRadians(rotationAngle), textLoc)); + anns.add(ann); + + Collection> bundles = stitch ? List.of(allCombined) : byParent.values(); + + for (List bundle : bundles) { + Location middle = findMiddleLoc(bundle); + + // don't draw the arrow all the way, leave a buffer on each end + double dist = LocationUtils.horzDistanceFast(middle, textLoc); + double az = LocationUtils.azimuthRad(textLoc, middle); + + double buffer; + if (dist > 40) + buffer = 10d; + else if (dist > 20) + buffer = 5; + else if (dist > 10) + buffer = 2; + else + buffer = 0.15*dist; + + Location start = LocationUtils.location(textLoc, az, buffer); + Location end = LocationUtils.location(textLoc, az, dist-buffer); + + arrows.add(LocationList.of(start, end)); + } + } + + private static Location findMiddleLoc(List sects) { + FaultTrace combTrace = new FaultTrace(null); + + for (FaultSection sect : sects) + combTrace.addAll(sect.getFaultTrace()); + + combTrace = FaultUtils.resampleTrace(combTrace, 2); + return combTrace.get(1); + } + +} diff --git a/src/main/java/scratch/kevin/prvi25/figures/DefModelSampleLineIntegralsPlot.java b/src/main/java/scratch/kevin/prvi25/figures/DefModelSampleLineIntegralsPlot.java index 7f411112..ff3e77dc 100644 --- a/src/main/java/scratch/kevin/prvi25/figures/DefModelSampleLineIntegralsPlot.java +++ b/src/main/java/scratch/kevin/prvi25/figures/DefModelSampleLineIntegralsPlot.java @@ -43,18 +43,22 @@ public class DefModelSampleLineIntegralsPlot { public static void main(String[] args) throws IOException { - File outputDir = new File("/tmp"); +// File outputDir = new File("/tmp"); + File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_ERF/prvi25-erf-paper/Figures/crustal_dm"); String prefix = "prvi_dm_sample_line_integrals"; LogicTree randTree = LogicTree.read(new File( "/home/kevin/OpenSHA/nshm23/batch_inversions/2024_10_24-prvi25_crustal_branches-dmSample5x/logic_tree.json")); PRVI25_CrustalFaultModels fm = PRVI25_CrustalFaultModels.PRVI_CRUSTAL_FM_V1p1; - VectorComponent[] components = VectorComponent.values(); +// VectorComponent[] components = VectorComponent.values(); + VectorComponent[] components = { VectorComponent.FULL_HORIZONTAL }; double[] fractiles = {0d, 0.025, 0.16, 0.84, 0.975d, 1d}; String fractileLabel = "Sample p[0,2.5,16,84,97.5,100]"; + boolean titles = false; + boolean plotOrig = true; CPT randColorCPT = GMT_CPT_Files.CATEGORICAL_TAB10.instance(); @@ -153,7 +157,7 @@ else if (fractiles[f] == 1) // Color distColor = Colors.tab_red; Color distColor = Colors.tab_grey; Color origColor = Colors.tab_green; - Color transColor = new Color(distColor.getRed(), distColor.getGreen(), distColor.getBlue(), 30); + Color transColor = new Color(distColor.getRed(), distColor.getGreen(), distColor.getBlue(), 60); PlotCurveCharacterstics minMaxChar = new PlotCurveCharacterstics(PlotLineType.SHADED_UNCERTAIN, 1f, transColor); funcs.add(bounds); @@ -171,6 +175,10 @@ else if (fractiles[f] == 1) funcs.add(lowerIntegrals[c]); chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, Color.BLACK)); + avgIntegral.setName("Average"); + funcs.add(avgIntegral); + chars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 3f, Color.BLACK)); + upperIntegrals[c].setName(null); funcs.add(upperIntegrals[c]); chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, Color.BLACK)); @@ -203,7 +211,7 @@ else if (fractiles[f] == 1) chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, upperColor)); } - PlotSpec plot = new PlotSpec(funcs, chars, "Deformation Model Line Integrals", "Longitude", + PlotSpec plot = new PlotSpec(funcs, chars, titles ? "Deformation Model Line Integrals" : " ", "Longitude", components[c].label+" Summed Rate (mm/yr)"); plot.setLegendInset(RectangleAnchor.TOP_RIGHT, 0.975, 0.975, 0.9, false); @@ -254,17 +262,23 @@ private static DiscretizedFunc getAvg(DiscretizedFunc min, DiscretizedFunc max) return ret; } - private static double minLat = 16; - private static double maxLat = 21; - private static double minLon = -71d; - private static double maxLon = -60.5d; +// private static double minLat = 16; +// private static double maxLat = 21; +// private static double minLon = -71d; +// private static double maxLon = -60.5d; + + private static double minLat = SlipRateFigures.CRUSTAL_FAULT_MAP_REG.getMinLat(); + private static double maxLat = SlipRateFigures.CRUSTAL_FAULT_MAP_REG.getMaxLat(); + private static double minLon = SlipRateFigures.CRUSTAL_FAULT_MAP_REG.getMinLon(); + private static double maxLon = SlipRateFigures.CRUSTAL_FAULT_MAP_REG.getMaxLon(); private static List calcLineIntegrals(FaultSystemLineIntegralCalculator calc, double delta) { List integrals = new ArrayList<>(); double minLon = DefModelSampleLineIntegralsPlot.minLon; if (delta == 1d) - minLon += 1d; +// minLon += 1d; + minLon = Math.ceil(minLon); for (double lon=minLon; (float)lon<=(float)maxLon; lon += delta) integrals.add(calc.calcLineIntegral(new Location(minLat, lon), new Location(maxLat, lon))); diff --git a/src/main/java/scratch/kevin/prvi25/figures/GridCubeRatePlot.java b/src/main/java/scratch/kevin/prvi25/figures/GridCubeRatePlot.java new file mode 100644 index 00000000..e271e84c --- /dev/null +++ b/src/main/java/scratch/kevin/prvi25/figures/GridCubeRatePlot.java @@ -0,0 +1,469 @@ +package scratch.kevin.prvi25.figures; + +import java.awt.Color; +import java.awt.Font; +import java.io.File; +import java.io.IOException; +import java.text.DecimalFormat; +import java.util.ArrayList; +import java.util.List; +import java.util.Map; +import java.util.function.Function; + +import org.apache.commons.math3.stat.StatUtils; +import org.jfree.chart.annotations.XYTextAnnotation; +import org.jfree.chart.ui.TextAnchor; +import org.jfree.data.Range; +import org.opensha.commons.data.function.DefaultXY_DataSet; +import org.opensha.commons.data.function.EvenlyDiscretizedFunc; +import org.opensha.commons.data.function.HistogramFunction; +import org.opensha.commons.data.function.XY_DataSet; +import org.opensha.commons.data.xyz.EvenlyDiscrXYZ_DataSet; +import org.opensha.commons.geo.CubedGriddedRegion; +import org.opensha.commons.geo.GriddedRegion; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.LocationUtils; +import org.opensha.commons.geo.Region; +import org.opensha.commons.geo.json.Feature; +import org.opensha.commons.geo.json.FeatureProperties; +import org.opensha.commons.geo.json.Geometry; +import org.opensha.commons.geo.json.Geometry.LineString; +import org.opensha.commons.geo.json.Geometry.GeometryCollection; +import org.opensha.commons.geo.json.Geometry.Polygon; +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.PlotUtils; +import org.opensha.commons.gui.plot.jfreechart.xyzPlot.XYZPlotSpec; +import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; +import org.opensha.commons.util.cpt.CPT; +import org.opensha.sha.earthquake.faultSysSolution.FaultSystemRupSet; +import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; +import org.opensha.sha.earthquake.faultSysSolution.util.FaultSysTools; +import org.opensha.sha.earthquake.rupForecastImpl.nshm23.gridded.NSHM23_FaultCubeAssociations; +import org.opensha.sha.earthquake.rupForecastImpl.nshm23.gridded.NSHM23_SingleRegionGridSourceProvider; +import org.opensha.sha.faultSurface.GeoJSONFaultSection; +import org.opensha.sha.magdist.GutenbergRichterMagFreqDist; +import org.opensha.sha.magdist.IncrementalMagFreqDist; +import org.opensha.sha.util.TectonicRegionType; + +import com.google.common.base.Preconditions; +import com.google.common.primitives.Doubles; + +import scratch.UCERF3.erf.ETAS.SeisDepthDistribution; + +public class GridCubeRatePlot { + + public static void main(String[] args) throws IOException { + double lonSpacingKM = LocationUtils.horzDistance(new Location(0d, 0d), new Location(0d, 1d)); + double gridSpacing = 0.1; + double cellSpacingKM = gridSpacing * lonSpacingKM; + double upDepth = 0d; + double lowDepth = 15d; +// int cellsBefore = 2; +// int cellsAfter = 2; +// double shadeInCellFract = 0.3; + int cellsBefore = 2; + int cellsAfter = 2; + double shadeInCellFract = 0.0; +// int cellsBefore = 1; +// int cellsAfter = 2; +// double shadeInCellFract = 0.0; + boolean vertical = false; + double dip, rake; + Location traceStart, traceEnd; + if (vertical) { + dip = 90d; + rake = 0d; + traceStart = new Location(0d, 0d); + traceEnd = new Location(1d, 0d); + } else { + // set the dip such that the association 2 cells away is zero + double assocWidth = NSHM23_SingleRegionGridSourceProvider.DEFAULT_MAX_FAULT_NUCL_DIST; + double fromSideOfCell = assocWidth % cellSpacingKM; + System.out.println("Calculated dist from site of cell: "+(float)fromSideOfCell); + Preconditions.checkState(fromSideOfCell*2 < cellSpacingKM); + double horzWidth = cellSpacingKM - 2*fromSideOfCell; + double height = lowDepth - upDepth; + dip = Math.toDegrees(Math.atan(height/horzWidth)); + double traceLon = -0.5*gridSpacing + fromSideOfCell/lonSpacingKM; + traceStart = new Location(0d, traceLon); + traceEnd = new Location(1d, traceLon); + rake = 0d; + } + + File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_ERF/prvi25-erf-paper/Figures/crustal_grid"); + + Region sectReg = new Region(new Location(traceStart.lat, traceStart.lon-1d), new Location(traceEnd.lat, traceEnd.lon+1d)); + GriddedRegion gridReg = new GriddedRegion(sectReg, gridSpacing, GriddedRegion.ANCHOR_0_0); + + GutenbergRichterMagFreqDist supraMFD = new GutenbergRichterMagFreqDist(1d, 0.1, 6.55, 7.55, 11); + EvenlyDiscretizedFunc gridRefMFD = FaultSysTools.initEmptyMFD(5.05, 7.55); + GutenbergRichterMagFreqDist gridMFD = new GutenbergRichterMagFreqDist(1d, 1d, gridRefMFD.getMinX(), gridRefMFD.getMaxX(), gridRefMFD.size()); + System.out.println("Have "+gridReg.getNodeCount()+" total grid nodes"); + double rateEach = 1e-2; + double sumRate = rateEach*gridReg.getNodeCount(); + gridMFD.scaleToCumRate(5.05, sumRate); + + FeatureProperties props = new FeatureProperties(); + props.set(GeoJSONFaultSection.FAULT_ID, 0); + props.set(GeoJSONFaultSection.FAULT_NAME, "Test Fault"); + props.set(GeoJSONFaultSection.DIP, dip); + props.set(GeoJSONFaultSection.RAKE, rake); + props.set(GeoJSONFaultSection.LOW_DEPTH, lowDepth); + props.set(GeoJSONFaultSection.UPPER_DEPTH, upDepth); + props.set(GeoJSONFaultSection.SLIP_RATE, 10d); + LineString geom = new LineString(traceStart, traceEnd); +// props.set(GeoJSONFaultSection.PROXY, true); +// Geometry geom = new LineString(traceStart, traceEnd); +// Polygon proxyZone = new Polygon(new Region(traceStart, new Location(traceEnd.lat, traceEnd.lon+gridSpacing))); +// geom = new GeometryCollection(geom, proxyZone); + Feature sectFeature = new Feature(geom, props); + GeoJSONFaultSection sect = GeoJSONFaultSection.fromFeature(sectFeature); + + Location centerLoc = new Location(0.5*(sectReg.getMinLat()+sectReg.getMaxLat()), + 0.5*(sectReg.getMinLon()+sectReg.getMaxLon())); + int centerGridIndex = gridReg.indexForLocation(centerLoc); + + List plotCells = new ArrayList<>(); + for (int i=0; i cellXRanges = new ArrayList<>(); + for (int cell : plotCells) { + Location loc = gridReg.getLocation(cell); + System.out.println("\t"+loc); + double cellStart = loc.lon-0.5*gridSpacing; + double cellEnd = loc.lon+0.5*gridSpacing; + plotMinLon = Math.min(plotMinLon, cellStart); + plotMaxLon = Math.max(plotMaxLon, cellEnd); + cellXRanges.add(new Range(cellStart*lonSpacingKM, cellEnd*lonSpacingKM)); + } + + Range xRange = new Range(plotMinLon*lonSpacingKM + shadeInCellFract*cellSpacingKM, plotMaxLon*lonSpacingKM - shadeInCellFract*cellSpacingKM); + + double[] rupMags = new double[supraMFD.size()]; + double[] rupAreas = new double[supraMFD.size()]; + double[] rupLengths = new double[supraMFD.size()]; + double[] rupRates = new double[supraMFD.size()]; + List> sectsForRup = new ArrayList<>(rupRates.length); + for (int i=0; i cellCubes = new ArrayList<>(); + double minCubeX = Double.POSITIVE_INFINITY; + double maxCubeX = Double.NEGATIVE_INFINITY; + double minCubeZ = Double.POSITIVE_INFINITY; + double maxCubeZ = Double.NEGATIVE_INFINITY; + for (int cell : plotCells) { + int[] cubes = cgr.getCubeIndicesForGridCell(cell); + cellCubes.add(cubes); + for (int cube : cubes) { + Location cubeLoc = cgr.getCubeLocationForIndex(cube); + double x = cubeLoc.lon*lonSpacingKM; + double z = cubeLoc.depth; + minCubeX = Math.min(minCubeX, x); + maxCubeX = Math.max(maxCubeX, x); + minCubeZ = Math.min(minCubeZ, z); + maxCubeZ = Math.max(maxCubeZ, z); + } + } + System.out.println("Cube xSpacing: "+(float)cubeSpacingX); + System.out.println("Cube zSpacing: "+(float)cubeSpacingZ); + System.out.println("Cube xRange: ["+(float)minCubeX+", "+(float)maxCubeX+"]"); + System.out.println("Cube zRange: ["+(float)minCubeZ+", "+(float)maxCubeZ+"]"); + int nx = (int)((maxCubeX - minCubeX)/cubeSpacingX) + 1; + int nz = (int)((maxCubeZ - minCubeZ)/cubeSpacingZ) + 1; + EvenlyDiscrXYZ_DataSet assocXYZ = new EvenlyDiscrXYZ_DataSet(nx, nz, minCubeX, minCubeZ, cubeSpacingX, cubeSpacingZ); + Preconditions.checkState((float)assocXYZ.getMaxX() == (float)maxCubeX, + "x mismatch; nx=%s, cubeMaxX=%s, xyzMaxX=%s", nx, maxCubeX, assocXYZ.getMaxX()); + Preconditions.checkState((float)assocXYZ.getMaxY() == (float)maxCubeZ, + "z mismatch; nz=%s, cubeMaxZ=%s, xyzMaxY=%s", nz, maxCubeZ, assocXYZ.getMaxY()); + List> xyzMappedCubes = new ArrayList<>(); + for (int i=0; i()); + for (int[] cubes : cellCubes) { + for (int cube : cubes) { + Location cubeLoc = cgr.getCubeLocationForIndex(cube); + double x = cubeLoc.lon*lonSpacingKM; + double z = cubeLoc.depth; + int xyzIndex = assocXYZ.indexOf(x, z); + Preconditions.checkArgument(xyzIndex >= 0); + xyzMappedCubes.get(xyzIndex).add(cube); + } + } + int minNumMapped = Integer.MAX_VALUE; + int maxNumMapped = 0; + for (List mapped : xyzMappedCubes) { + minNumMapped = Integer.min(minNumMapped, mapped.size()); + maxNumMapped = Integer.max(maxNumMapped, mapped.size()); + } + System.out.println("Cube to xyz mapping range: ["+minNumMapped+", "+maxNumMapped+"]"); + + XY_DataSet sectXY = new DefaultXY_DataSet(); + Location sectTopLoc = sect.getFaultTrace().first(); + Location sectBotLoc = sect.getFaultSurface(1d).getEvenlyDiscritizedLowerEdge().first(); + sectXY.set(sectTopLoc.lon*lonSpacingKM, sectTopLoc.depth); + sectXY.set(sectBotLoc.lon*lonSpacingKM, sectBotLoc.depth); + + CPT assocCPT = GMT_CPT_Files.SEQUENTIAL_LAJOLLA_UNIFORM.instance().rescale(0d, 1d); + assocCPT.setPreferredTickInterval(0.1); + calcXYZValues(assocXYZ, xyzMappedCubes, new Function() { + + @Override + public Double apply(Integer cubeIndex) { + double[] weights = cubeAssoc.getOrigSectDistWeightsAtCube(cubeIndex); + if (weights == null) + return 0d; + return StatUtils.sum(weights); + } + }, false); // average +// List cellAssocs = aggregateToCells(assocXYZ, cellXRanges, false); + List cellAssocs = new ArrayList<>(); + for (int cell : plotCells) + cellAssocs.add(cubeAssoc.getNodeFraction(cell)); + System.out.println("Cell associations: "+cellAssocs); + + DecimalFormat fractDF = new DecimalFormat("0.00"); + + XYZPlotSpec assocPlot = buildPlot(assocXYZ, cellXRanges, cellAssocs, new FormatFunc(fractDF), + sectXY, assocCPT, "Fractional Association", " "); + + HeadlessGraphPanel gp = PlotUtils.initHeadless(); + + Range yRange = new Range(CELL_TOP_Z, assocXYZ.getMaxY()+0.5*cubeSpacingZ); + gp.drawGraphPanel(assocPlot, false, false, xRange, yRange); + + int width = plotCells.size() > 4 ? 1000 : 800; + + PlotUtils.writePlots(outputDir, "cube_assoc", gp, width, false, true, true, false); + + EvenlyDiscrXYZ_DataSet totalNuclXYZ = assocXYZ.copy(); + EvenlyDiscrXYZ_DataSet m7NuclXYZ = assocXYZ.copy(); + calcXYZValues(totalNuclXYZ, xyzMappedCubes, new Function() { + + @Override + public Double apply(Integer cubeIndex) { +// return gridProv.getTotalMFD_ForCube(cubeIndex).calcSumOfY_Vals(); + return gridProv.getGriddedSeisMFD_ForCube(cubeIndex).calcSumOfY_Vals(); + } + }, true); // sum them + calcXYZValues(m7NuclXYZ, xyzMappedCubes, new Function() { + + @Override + public Double apply(Integer cubeIndex) { +// return gridProv.getTotalMFD_ForCube(cubeIndex).getCumRate(gridRefMFD.getClosestXIndex(7.01)); + return gridProv.getGriddedSeisMFD_ForCube(cubeIndex).getCumRate(gridRefMFD.getClosestXIndex(7.01)); + } + }, true); // sum them + List totCellRates = new ArrayList<>(); + List m7CellRates = new ArrayList<>(); + totalNuclXYZ.log10(); + m7NuclXYZ.log10(); + for (int cell : plotCells) { + IncrementalMagFreqDist cellMFD = gridProv.getMFD(cell); + totCellRates.add(cellMFD.calcSumOfY_Vals()); + m7CellRates.add(cellMFD.getCumRate(cellMFD.getClosestXIndex(7.01))); + } + System.out.println("Cell total nucl rates: "+totCellRates); + System.out.println("Cell M>7 nucl rates: "+m7CellRates); + for (int i=0; i rateFormat = new Function() { + + private DecimalFormat eDF = new DecimalFormat("0.000E0"); + + @Override + public String apply(Double t) { + // back to linear + double linear = Math.pow(10, t); + return eDF.format(linear).replace("E", "e"); + } + }; + + XYZPlotSpec cellRatePlot = buildPlot(totalNuclXYZ, cellXRanges, totCellRates, rateFormat, + sectXY, nuclRateCPT, "Gridded M>5 Nucleation Rate", " "); + + gp.drawGraphPanel(cellRatePlot, false, false, xRange, yRange); + + PlotUtils.writePlots(outputDir, "nucl_rate_tot", gp, width, false, true, true, false); + + XYZPlotSpec cellRateM7Plot = buildPlot(m7NuclXYZ, cellXRanges, m7CellRates, rateFormat, + sectXY, nuclRateCPT, "Gridded M>7 Nucleation Rate", " "); + + gp.drawGraphPanel(cellRateM7Plot, false, false, xRange, yRange); + + PlotUtils.writePlots(outputDir, "nucl_rate_m7", gp, width, false, true, true, false); + } + + private static void calcXYZValues(EvenlyDiscrXYZ_DataSet xyz, List> xyzMappedCubes, + Function cubeValFunc, boolean sum) { + for (int i=0; i mappedCubes = xyzMappedCubes.get(i); + Preconditions.checkState(!mappedCubes.isEmpty()); + double sumVals = 0d; + for (int cube : mappedCubes) + sumVals += cubeValFunc.apply(cube); + if (sum) + xyz.set(i, sumVals); + else + // take average + xyz.set(i, sumVals/mappedCubes.size()); + } + } + + private static final double CELL_TOP_Z = -1.8d; + private static final double CELL_LINE_SHIFT = 0.05; + private static final double CELL_LINE_TOP_Z = CELL_TOP_Z+CELL_LINE_SHIFT; + private static final double CELL_LINE_BOT_Z = -CELL_LINE_SHIFT; + + private static List aggregateToCells(EvenlyDiscrXYZ_DataSet xyz, List cellXRanges, boolean sum) { + double[] cellSums = new double[cellXRanges.size()]; + int[] cellCounts = new int[cellXRanges.size()]; + for (int i=0; i= 0, "Cube with x=%s not contained in cellXRanges=%s", x, cellXRanges); + cellSums[cellIndex] += val; + cellCounts[cellIndex]++; + } + if (!sum) { + // take average + List ret = new ArrayList<>(); + for (int i=0; i 0); + ret.add(cellSums[i]/cellCounts[i]); + } + return ret; + } + return Doubles.asList(cellSums); + } + + private static class FormatFunc implements Function { + + private DecimalFormat df; + + public FormatFunc(DecimalFormat df) { + this.df = df; + } + + @Override + public String apply(Double t) { + return df.format(t); + } + + } + + private static XYZPlotSpec buildPlot(EvenlyDiscrXYZ_DataSet xyz, List cellXRanges, List cellValues, + Function cellDF, XY_DataSet sectXY, CPT cpt, String zLabel, String title) { + List funcs = new ArrayList<>(); + List chars = new ArrayList<>(); + + List cellAnns = cellDF == null ? null : new ArrayList<>(); + Font annFont = new Font(Font.SANS_SERIF, Font.BOLD, cellValues.size() > 3 ? 16 : 20); + + for (int i=0; i funcs = new ArrayList<>(); + List chars = new ArrayList<>(); + + if (obs != null) { + funcs.add(obs); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, obsColor)); + + IncrementalMagFreqDist obsLower = obs.getLower(); + obsLower.setName(obs.getBoundName()); + funcs.add(obsLower); + chars.add(new PlotCurveCharacterstics(PlotLineType.DOTTED, 2f, obsColor)); + + IncrementalMagFreqDist obsUpper = obs.getUpper(); + obsUpper.setName(null); + funcs.add(obsUpper); + chars.add(new PlotCurveCharacterstics(PlotLineType.DOTTED, 2f, obsColor)); + } + + onFaultMean.setName("On-Fault"); + funcs.add(onFaultMean); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, onFaultColor)); + + if (onFaultDists != null) { + for (int i=0; i sectIDs = new ArrayList<>(); + if (parents != null && parents.length > 0) { + for (FaultSection sect : sol.getRupSet().getFaultSectionDataList()) + if (Ints.contains(parents, sect.getParentSectionId())) + sectIDs.add(sect.getSectionId()); + } else { + for (int i=0; i nameRemappings = new HashMap<>(); + nameRemappings.put(NSHM23_ScalingRelationships.WIDTH_LIMITED.getShortName(), "Wdth-Lmtd"); + nameRemappings.put(NSHM23_ScalingRelationships.WIDTH_LIMITED_CSD.getShortName(), "W-L, CSD"); + nameRemappings.put(NSHM23_ScalingRelationships.LOGA_C4p2_SQRT_LEN.getShortName(), "Sqrt-Len"); + nameRemappings.put(PRVI25_CrustalRandomlySampledDeformationModelLevel.NAME, "Geologic Deformation Model Sample"); + + File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_ERF/prvi25-erf-paper/Figures/logic_trees"); + + List> trees = new ArrayList<>(); + List prefixes = new ArrayList<>(); + + trees.add(LogicTree.read(new File("/data/kevin/nshm23/batch_inversions/" + + "2024_11_19-prvi25_crustal_branches-dmSample5x/logic_tree.json"))); +// trees.set(trees.size()-1, reorderTree(trees.get(trees.size()-1), PRVI25_CrustalRandomlySampledDeformationModelLevel.NAME, 0)); + prefixes.add("crustal_inversion"); + trees.add(new PRVI25_InvConfigFactory().getGridSourceTree(trees.get(trees.size()-1))); + prefixes.add("crustal_gridded"); + + trees.add(LogicTree.read(new File("/data/kevin/nshm23/batch_inversions/" + + "2024_11_19-prvi25_subduction_branches/logic_tree.json"))); + prefixes.add("subduction_inversion"); + trees.add(new PRVI25_InvConfigFactory().getGridSourceTree(trees.get(trees.size()-1))); + prefixes.add("subduction_gridded"); + + for (int i=0; i tree = trees.get(i); + String prefix = prefixes.get(i); + System.out.println("Plotting "+prefix); + + LogicTreeFigureWriter panel = new LogicTreeFigureWriter(tree, false, nameRemappings); + panel.write(outputDir, prefix, true, true); + + System.out.println(); + } + } + + private static LogicTree reorderTree(LogicTree tree, String levelName, int targetIndex) { + int foundIndex = -1; + List> newLevels = new ArrayList<>(); + for (int l=0; l level = tree.getLevels().get(l); + if (level.getName().equals(levelName)) + foundIndex = l; + newLevels.add(level); + } + Preconditions.checkState(foundIndex >= 0, "Didn't find level: %s", levelName); + newLevels.add(targetIndex, newLevels.remove(foundIndex)); + List> newBranches = new ArrayList<>(tree.size()); + for (LogicTreeBranch branch : tree) { + LogicTreeBranch newBranch = new LogicTreeBranch<>(newLevels); + for (LogicTreeNode node : branch) + newBranch.setValue(node); + newBranches.add(newBranch); + } + + return LogicTree.fromExisting(newLevels, newBranches); + } + + private static final int levelGap = 10; + private static final int levelHGap = 10; + private static final int choiceHGap = 5; + private static final int lineGap = 5; + private static final int minWidthPerNode = 100; + private static final Font levelFont = new Font(Font.SANS_SERIF, Font.BOLD, 24); + private static final Font choiceFont = new Font(Font.SANS_SERIF, Font.PLAIN, 18); + private static final Font weightFont = new Font(Font.SANS_SERIF, Font.ITALIC, 18); + private static final DecimalFormat weightDF = new DecimalFormat("0.###"); + private int levelFontHeight; + private int choiceFontHeight; + private int choiceWidth; + private int lineHeight; + private int levelHeight; + private int width; + private int height; + private List> includedLevels; + private List> includedLevelChoices; + + private LogicTree tree; + private Map nameRemappings; + private double totalWeight; + private Map choiceWeights; + private int maxNodes; + + public LogicTreeFigureWriter(LogicTree tree, boolean includeSingleChoice, Map nameRemappings) { + this.tree = tree; + if (nameRemappings == null) + nameRemappings = Map.of(); + this.nameRemappings = nameRemappings; + + maxNodes = 1; + int maxChoiceFontWidth = 0; + int maxLevelFontWidth = 0; + includedLevels = new ArrayList<>(); + includedLevelChoices = new ArrayList<>(); + totalWeight = 0d; + choiceWeights = new HashMap<>(); + for (int i=0; i level : tree.getLevels()) { + // figure out how many choices I have + List uniqueChoices = new ArrayList<>(); + HashSet uniqueChoicesSet = new HashSet<>(); + for (LogicTreeBranch branch : tree) { + LogicTreeNode node = branch.requireValue(level.getType()); + double weight = tree.getBranchWeight(branch); + if (!uniqueChoicesSet.contains(node)) { + uniqueChoices.add(node); + uniqueChoicesSet.add(node); + choiceWeights.put(node, weight); + } else { + choiceWeights.put(node, choiceWeights.get(node)+weight); + } + } + System.out.println(level.getName()+" has "+uniqueChoices.size()+" unique nodes"); + if (includeSingleChoice || uniqueChoices.size() > 1) { + includedLevels.add(level); + includedLevelChoices.add(uniqueChoices); + maxLevelFontWidth = Integer.max(maxLevelFontWidth, new JPanel().getFontMetrics(levelFont).stringWidth(remapped(level.getName()))); + if (!(level instanceof RandomlySampledLevel)) { + maxNodes = Integer.max(maxNodes, uniqueChoices.size()); + for (LogicTreeNode node : uniqueChoices) { + String name = remapped(node.getShortName()); + FontMetrics metrics = new JPanel().getFontMetrics(choiceFont); + maxChoiceFontWidth = Integer.max(maxChoiceFontWidth, metrics.stringWidth(name)); + } + } + } + } + + System.out.println("We have "+includedLevelChoices.size()+" levels and at most "+maxNodes+" nodes"); + + FontMetrics metrics = new JPanel().getFontMetrics(levelFont); + levelFontHeight = metrics.getHeight(); + System.out.println("Level font height in pixels: " + levelFontHeight); + metrics = new JPanel().getFontMetrics(choiceFont); + choiceFontHeight = metrics.getHeight(); + System.out.println("Choice font height in pixels: " + choiceFontHeight); + lineHeight = (int)(1.5d*levelFontHeight + 0.5); + levelHeight = levelFontHeight + lineGap*2 + lineHeight + choiceFontHeight + choiceFontHeight; + + choiceWidth = Integer.max(minWidthPerNode, maxChoiceFontWidth+choiceHGap); + width = Integer.max((int)(choiceWidth*maxNodes + 0.5*choiceWidth + 0.5), + maxLevelFontWidth+2*levelHGap); + height = levelHeight * includedLevelChoices.size() + levelGap * includedLevels.size(); + + System.out.println("Calculated dimensions: "+width+"x"+height); + + setSize(width, height); + setBackground(Color.WHITE); + } + + private String remapped(String name) { + if (nameRemappings.containsKey(name)) + return nameRemappings.get(name); + return name; + } + + @Override + protected void paintComponent(Graphics g) { + super.paintComponent(g); + + int y = (int)(0.5*levelGap + 0.5); + int middleX = (int)(0.5*width + 0.5); + + Graphics2D g2d = (Graphics2D) g; + float lineWidth = 3f; + BasicStroke lineStroke = new BasicStroke(lineWidth); + BasicStroke randLineStroke = new BasicStroke(lineWidth, BasicStroke.CAP_BUTT, + BasicStroke.JOIN_BEVEL,0,new float[] {Float.min(6, Float.max(lineWidth*0.7f, 3))},0); + g2d.setStroke(new BasicStroke(3)); + + for (int l=0; l level = includedLevels.get(l); + String name = remapped(level.getName()); + g2d.setFont(levelFont); + drawText(g2d, name, middleX, y, RectangleAnchor.TOP); + + y += levelFontHeight; + y += lineGap; + List nodes = includedLevelChoices.get(l); + int botLineY = y+lineHeight; + int topChoiceY = botLineY + lineGap; + int topWeightY = topChoiceY + choiceFontHeight; + if (level instanceof RandomlySampledLevel && nodes.size() > maxNodes && includedLevels.size() > 1) { + // figure out how many branches we have without this + HashSet uniquesWithout = new HashSet<>(); + for (LogicTreeBranch branch : tree) { + String str = ""; + for (int l1=0; l1 maxNodes) + prefNumLines = maxNodes; + int myNodesWidth = choiceWidth*prefNumLines; + int sideBuffer = (int)(0.5*(width - myNodesWidth)); + int leftX = sideBuffer; + g2d.setStroke(randLineStroke); + for (int i=0; i funcs = new ArrayList<>(); + List chars = new ArrayList<>(); + + Double m1 = null; + Double mMax = null; + + EvenlyDiscretizedFunc refMFD = FaultSysTools.initEmptyMFD(3.05, 10.95); + + Color[] colors = { + Colors.tab_blue, + Colors.tab_orange, + Colors.tab_green + }; + + for (RateType type : RateType.values()) { + List rates = PRVI25_RegionalSeismicity.loadRates(PRVI25_SeismicityRegions.CRUSTAL, type); + + RateRecord meanRec = SeismicityRateFileLoader.locateMean(rates); + if (m1 == null) + m1 = meanRec.M1; + if (mMax == null && meanRec instanceof PureGR) + mMax = ((PureGR)meanRec).Mmax; + + if (funcs.isEmpty()) { + EvenlyDiscretizedFunc mfd; + if (incremental) + mfd = SeismicityRateFileLoader.buildIncrementalMFD(meanRec, refMFD, refMFD.getMaxX()); + else + mfd = cmlMFD(meanRec, refMFD); + mfd.setName("Mean"); + funcs.add(mfd); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 4f, Color.BLACK)); + } + + Color color = colors[type.ordinal()]; + RateRecord low = SeismicityRateFileLoader.locateQuantile(rates, 0.025); + RateRecord high = SeismicityRateFileLoader.locateQuantile(rates, 0.975); + + EvenlyDiscretizedFunc lowMFD; + EvenlyDiscretizedFunc highMFD; + if (incremental) { + lowMFD = SeismicityRateFileLoader.buildIncrementalMFD(low, refMFD, refMFD.getMaxX()); + highMFD = SeismicityRateFileLoader.buildIncrementalMFD(high, refMFD, refMFD.getMaxX()); + } else { + lowMFD = cmlMFD(low, refMFD); + highMFD = cmlMFD(high, refMFD); + } + + lowMFD.setName(type.toString()); + funcs.add(lowMFD); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, color)); + highMFD.setName(null); + funcs.add(highMFD); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, color)); + } + + for (XY_DataSet func : funcs) { + if (func.getName() != null && func.getName().contains("M1")) + func.setName(func.getName().replace("M1", "M₁")); + if (func.getName() != null && func.getName().contains("Mmax")) + func.setName(func.getName().replace("Mmax", "Mₘₐₓ")); + } + + Range xRange = new Range(4d, 8d); + Range yRange = incremental ? new Range(1e-5, 1e1) : new Range(1e-4, 1e2); + + List anns = new ArrayList<>(); + Font annFont = new Font(Font.SANS_SERIF, Font.PLAIN, 22); + + DefaultXY_DataSet m1Line = new DefaultXY_DataSet(); + m1Line.set(m1, yRange.getLowerBound()); + m1Line.set(m1, yRange.getUpperBound()); + funcs.add(m1Line); + chars.add(new PlotCurveCharacterstics(PlotLineType.DOTTED, 2f, Color.DARK_GRAY)); + + DefaultXY_DataSet mMaxLine = new DefaultXY_DataSet(); + mMaxLine.set(mMax, yRange.getLowerBound()); + mMaxLine.set(mMax, yRange.getUpperBound()); + funcs.add(mMaxLine); + chars.add(new PlotCurveCharacterstics(PlotLineType.DOTTED, 2f, Color.DARK_GRAY)); + + DefaultXY_DataSet mc2Line = new DefaultXY_DataSet(); + mc2Line.set(6d, yRange.getLowerBound()); + mc2Line.set(6d, yRange.getUpperBound()); + funcs.add(mc2Line); + chars.add(new PlotCurveCharacterstics(PlotLineType.DOTTED, 2f, Color.DARK_GRAY)); + + XYTextAnnotation m1Ann = new XYTextAnnotation(" M₁=5", m1, yRange.getUpperBound()); + m1Ann.setFont(annFont); + m1Ann.setTextAnchor(TextAnchor.TOP_LEFT); + anns.add(m1Ann); + + XYTextAnnotation mMaxAnn = new XYTextAnnotation("Mₘₐₓ="+mMax.floatValue()+" ", mMax, yRange.getLowerBound()); + mMaxAnn.setFont(annFont); + mMaxAnn.setTextAnchor(TextAnchor.BOTTOM_RIGHT); + anns.add(mMaxAnn); + + XYTextAnnotation mc1Ann = new XYTextAnnotation(" Mc₁₉₇₃=5", m1, yRange.getLowerBound()); + mc1Ann.setFont(annFont); + mc1Ann.setTextAnchor(TextAnchor.BOTTOM_LEFT); + anns.add(mc1Ann); + + XYTextAnnotation mc2Ann = new XYTextAnnotation(" Mc₁₉₀₀=6", 6d, yRange.getLowerBound()); + mc2Ann.setFont(annFont); + mc2Ann.setTextAnchor(TextAnchor.BOTTOM_LEFT); + anns.add(mc2Ann); + + PlotSpec plot = new PlotSpec(funcs, chars, " ", "Magnitude", incremental ? "Incremental Rate (1/yr)" : "Cumulative Rate (1/yr)"); + plot.setLegendInset(true); + plot.setPlotAnnotations(anns); + + HeadlessGraphPanel gp = PlotUtils.initHeadless(); + + gp.drawGraphPanel(plot, false, true, xRange, yRange); + + PlotUtils.writePlots(outputDir, prefix, gp, 700, 650, true, true, false); + } + + private static EvenlyDiscretizedFunc cmlMFD(RateRecord record, EvenlyDiscretizedFunc refMFD) { + if (record.type == RateType.EXACT) + return ((Exact)record).cumulativeDist; + Preconditions.checkState(record instanceof PureGR); + PureGR grRec = (PureGR)record; + GutenbergRichterMagFreqDist grMFD = new GutenbergRichterMagFreqDist( + grRec.b, 1d, refMFD.getMinX(), refMFD.getMaxX(), refMFD.size()); + grMFD.scaleToCumRate(grMFD.getX(grMFD.getClosestXIndex(grRec.M1+0.01)), grRec.rateAboveM1); + return grMFD.getCumRateDistWithOffset(); + } + +} diff --git a/src/main/java/scratch/kevin/prvi25/figures/PRVI_SubductionSubSectPlots.java b/src/main/java/scratch/kevin/prvi25/figures/PRVI_SubductionSubSectPlots.java index 92a8d564..988a7b76 100644 --- a/src/main/java/scratch/kevin/prvi25/figures/PRVI_SubductionSubSectPlots.java +++ b/src/main/java/scratch/kevin/prvi25/figures/PRVI_SubductionSubSectPlots.java @@ -16,6 +16,7 @@ import org.opensha.commons.geo.LocationList; import org.opensha.commons.geo.Region; import org.opensha.commons.gui.plot.GeographicMapMaker; +import org.opensha.commons.gui.plot.PlotUtils; import org.opensha.commons.logicTree.LogicTreeBranch; import org.opensha.commons.logicTree.LogicTreeNode; import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; @@ -48,7 +49,8 @@ public static void main(String[] args) throws IOException { branch.setValue(scale); PRVI25_InvConfigFactory factory = new PRVI25_InvConfigFactory(); - File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_Subduction/figures/fault_model"); +// File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_Subduction/figures/fault_model"); + File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_ERF/prvi25-erf-paper/Figures/sub_fm"); Preconditions.checkState(outputDir.exists() || outputDir.mkdir()); Font font = new Font(Font.SANS_SERIF, Font.BOLD, 24); @@ -58,6 +60,10 @@ public static void main(String[] args) throws IOException { CPT slipCPT = GMT_CPT_Files.SEQUENTIAL_BATLOW_UNIFORM.instance().rescale(0d, 5d); CPT slipUncertCPT = GMT_CPT_Files.SEQUENTIAL_BATLOW_UNIFORM.instance().rescale(0d, 2d); CPT rakeCPT = GMT_CPT_Files.SEQUENTIAL_NAVIA_UNIFORM.instance().rescale(0d, 90d); + + PlotUtils.writeScaleLegendOnly(outputDir, "slip_cpt", + GeographicMapMaker.buildCPTLegend(slipCPT, "Slip Rate (mm/yr)"), + GeographicMapMaker.PLOT_WIDTH_DEFAULT, true, true); for (PRVI25_SubductionFaultModels fm : PRVI25_SubductionFaultModels.values()) { if (fm.getNodeWeight(branch) == 0d) continue; @@ -71,6 +77,7 @@ public static void main(String[] args) throws IOException { List sects = rupSet.getFaultSectionDataList(); GeographicMapMaker mapMaker = new GeographicMapMaker(plotReg); + mapMaker.setWriteGeoJSON(false); mapMaker.setFaultSections(sects); mapMaker.setFillSurfaces(true); @@ -125,24 +132,32 @@ public static void main(String[] args) throws IOException { rakes.add(sect.getAveRake()); } - mapMaker.plotSectScalars(slips, slipCPT, dm.getShortName()+" Slip Rate (mm/yr)"); - mapMaker.plot(outputDir, "subduction_slip_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), fmName+", "+dmName); - - mapMaker.plotSectScalars(slipUncerts, slipUncertCPT, dm.getShortName()+" Slip Rate Uncertainty (mm/yr)"); - mapMaker.plot(outputDir, "subduction_slip_uncert_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), fmName+", "+dmName); +// mapMaker.plotSectScalars(slipUncerts, slipUncertCPT, dm.getShortName()+" Slip Rate Uncertainty (mm/yr)"); +// mapMaker.plot(outputDir, "subduction_slip_uncert_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), fmName+", "+dmName); - mapMaker.plotSectScalars(rakes, rakeCPT, dm.getShortName()+" Rake"); // draw rake lines int rakeMod = 3; List rakeArrows = new ArrayList<>(); for (FaultSection sect : sects) if (sect.getSectionId() % rakeMod == 0) rakeArrows.addAll(buildRakeArrows(sect, slipCPT.getMaxValue())); + mapMaker.plotArrows(rakeArrows, 20d, Color.BLACK, 2f); + mapMaker.setFillArrowheads(true); mapMaker.plotLines(rakeArrows, Color.BLACK, 2f); - mapMaker.plot(outputDir, "subduction_rake_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), fmName+", "+dmName); - mapMaker.clearLines(); + mapMaker.plotSectScalars(slips, slipCPT, dm.getShortName()+" Slip Rate (mm/yr)"); + mapMaker.plot(outputDir, "subduction_slip_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), fmName+", "+dmName); + + mapMaker.plotSectScalars(slips, slipCPT, null); + mapMaker.plot(outputDir, "subduction_slip_"+fm.getFilePrefix()+"_"+dm.getFilePrefix()+"_no_cpt", fmName+", "+dmName); + mapMaker.plotSectScalars(rakes, rakeCPT, dm.getShortName()+" Rake"); + String prefix = "subduction_rake_"+fm.getFilePrefix()+"_"+dm.getFilePrefix();; + System.out.println("Plotting "+prefix); + mapMaker.plot(outputDir, prefix, fmName+", "+dmName); + + mapMaker.clearLines(); + mapMaker.clearArrows(); } } } @@ -173,9 +188,9 @@ public static List buildRakeArrows(FaultSection sect, Location ref Location lineEnd = FaultSystemLineIntegralCalculator.locConstPlotDist(center, refLoc, az, len); rakeArrows.add(LocationList.of(center, lineEnd)); - Location arrowStart = FaultSystemLineIntegralCalculator.locConstPlotDist(lineEnd, center, az - 135, 0.2*len); - Location arrowEnd = FaultSystemLineIntegralCalculator.locConstPlotDist(lineEnd, center, az + 135, 0.2*len); - rakeArrows.add(LocationList.of(arrowStart, lineEnd, arrowEnd)); +// Location arrowStart = FaultSystemLineIntegralCalculator.locConstPlotDist(lineEnd, center, az - 135, 0.2*len); +// Location arrowEnd = FaultSystemLineIntegralCalculator.locConstPlotDist(lineEnd, center, az + 135, 0.2*len); +// rakeArrows.add(LocationList.of(arrowStart, lineEnd, arrowEnd)); return rakeArrows; } diff --git a/src/main/java/scratch/kevin/prvi25/figures/ProxyFaultRepresentationFigures.java b/src/main/java/scratch/kevin/prvi25/figures/ProxyFaultRepresentationFigures.java index 8bbb565a..7939e7b8 100644 --- a/src/main/java/scratch/kevin/prvi25/figures/ProxyFaultRepresentationFigures.java +++ b/src/main/java/scratch/kevin/prvi25/figures/ProxyFaultRepresentationFigures.java @@ -24,9 +24,9 @@ import org.opensha.commons.gui.plot.PlotCurveCharacterstics; import org.opensha.commons.gui.plot.PlotSymbol; import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; +import org.opensha.commons.util.DataUtils.MinMaxAveTracker; import org.opensha.commons.util.cpt.CPT; import org.opensha.sha.earthquake.faultSysSolution.FaultSystemRupSet; -import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; import org.opensha.sha.earthquake.faultSysSolution.modules.ProxyFaultSectionInstances; import org.opensha.sha.earthquake.faultSysSolution.ruptures.util.UniqueRupture; import org.opensha.sha.earthquake.faultSysSolution.util.FaultSectionUtils; @@ -57,192 +57,247 @@ public static void main(String[] args) throws IOException { FaultSystemRupSet rupSet = new PRVI25_InvConfigFactory().buildRuptureSet( PRVI25_LogicTreeBranch.DEFAULT_CRUSTAL_ON_FAULT, FaultSysTools.defaultNumThreads()); - rupSet.addModule(ProxyFaultSectionInstances.build(rupSet, 5, 5d)); + ProxyFaultSectionInstances proxyModule = ProxyFaultSectionInstances.build(rupSet, 5, 5d); + rupSet.addModule(proxyModule); -// int parentID = FaultSectionUtils.findParentSectionID(rupSet.getFaultSectionDataList(), "Anegada", "SW"); -// Region region = new Region(new Location(17.25, -66.25), new Location(18.5, -64)); -// int parentID = FaultSectionUtils.findParentSectionID(rupSet.getFaultSectionDataList(), "SW", "Puerto", "Rico"); -// Region region = new Region(new Location(17.7, -67.05), new Location(18.1, -66.6)); - int parentID = FaultSectionUtils.findParentSectionID(rupSet.getFaultSectionDataList(), "Cerro", "Goden"); - Region region = new Region(new Location(18.08, -67.4), new Location(18.4, -66.95)); + File outputDir = new File("/home/kevin/Documents/papers/2024_PRVI_ERF/prvi25-erf-paper/Figures/proxy_faults"); - String parentName = null; + boolean titles = false; - CPT colors = GMT_CPT_Files.CATEGORICAL_BATLOW_UNIFORM.instance(); +// int singleParent = FaultSectionUtils.findParentSectionID(rupSet.getFaultSectionDataList(), "Anegada", "SW"); + int singleParent = -1; - List proxySects = new ArrayList<>(); - for (FaultSection sect : rupSet.getFaultSectionDataList()) { - if (sect.getParentSectionId() == parentID) { - proxySects.add(sect); - if (parentName == null) - parentName = sect.getParentSectionName(); - } - } + double lengthLimit = 75d; - double lengthLimit = PRVI25_InvConfigFactory.MAX_PROXY_FAULT_RUP_LEN_DEFAULT*1e3; // km -> m +// RupType[] types = RupType.values(); + RupType[] types = { + RupType.PROXY, + RupType.SPLIT_PROXIES + }; - int rupID = -1; - double maxRupMag = 0d; - List allRups = new ArrayList<>(rupSet.getRupturesForParentSection(parentID)); - Collections.reverse(allRups); - for (int rupIndex : allRups) { - double mag = rupSet.getMagForRup(rupIndex); - if ((float)mag < (float)maxRupMag) + CPT colors = GMT_CPT_Files.CATEGORICAL_BATLOW_UNIFORM.instance(); + + Map parentProxyRegTracks = new HashMap<>(); + Map parentProxyNames = new HashMap<>(); + for (FaultSection sect : rupSet.getFaultSectionDataList()) { + Region poly = sect.getZonePolygon(); + if (!sect.isProxyFault() || poly == null) continue; - if (rupSet.getLengthForRup(rupIndex) > lengthLimit) + int parentID = sect.getParentSectionId(); + if (singleParent >= 0 && parentID != singleParent) continue; - List rupSects = rupSet.getFaultSectionDataForRupture(rupIndex); -// if (rupSects.size() != proxySects.size()) -// continue; - boolean allMatch = true; - for (FaultSection sect : rupSects) { - if (sect.getParentSectionId() != parentID) { - allMatch = false; - break; - } - } - if (allMatch) { - rupID = rupIndex; - maxRupMag = mag; - } - } - Preconditions.checkState(rupID >= 0); - - System.out.println(parentName); - System.out.println("Using an M"+(float)maxRupMag+" rupture, len="+(float)rupSet.getLengthForRup(rupID)*1e-3); - - File outputDiur = new File("/tmp/"); - - UniqueRupture rupUnique = UniqueRupture.forIDs(rupSet.getSectionsIndicesForRup(rupID)); - - LocationList gridLocs = new LocationList(); - for (FaultSection sect : proxySects) { - if (rupUnique.contains(sect.getSectionId())) { - GriddedRegion grid = new GriddedRegion(sect.getZonePolygon(), 0.05d, GriddedRegion.ANCHOR_0_0); - gridLocs.addAll(grid.getNodeList()); + + MinMaxAveTracker[] parentTracks = parentProxyRegTracks.get(parentID); + if (parentTracks == null) { + parentTracks = new MinMaxAveTracker[2]; + parentTracks[0] = new MinMaxAveTracker(); + parentTracks[1] = new MinMaxAveTracker(); + parentProxyRegTracks.put(parentID, parentTracks); + parentProxyNames.put(parentID, sect.getParentSectionName()); } + parentTracks[0].addValue(poly.getMinLat()); + parentTracks[0].addValue(poly.getMaxLat()); + parentTracks[1].addValue(poly.getMinLon()); + parentTracks[1].addValue(poly.getMaxLon()); +// for (FaultSection proxySect : proxyModule.getProxySects()) { +// if (proxySect.getParentSectionId() == sect.getSectionId()) { +// for (Location loc : proxySect.getFaultSurface(1d).getPerimeter()) { +// parentTracks[0].addValue(loc.lat); +// parentTracks[1].addValue(loc.lon); +// } +// } +// } } - ProxyFaultSectionInstances proxyModule = rupSet.requireModule(ProxyFaultSectionInstances.class); - - RuptureSurface rup = rupSet.getSurfaceForRupture(rupID, 1d); - - for (RupType type : RupType.values()) { - GeographicMapMaker mapMaker = new GeographicMapMaker(region); + for (int parentID : parentProxyNames.keySet()) { + String parentName = parentProxyNames.get(parentID); + MinMaxAveTracker[] regTracks = parentProxyRegTracks.get(parentID); + Location lowerLeft = new Location(regTracks[0].getMin(), regTracks[1].getMin()); + Location upperRight = new Location(regTracks[0].getMax(), regTracks[1].getMax()); - if (type == RupType.PROXY) { - mapMaker.setFaultSections(proxySects); - Color color = color(colors, 0); - List sectColors = new ArrayList<>(); - for (int i=0; i chars = new ArrayList<>(); - for (int i=0; i combFaultSects = new ArrayList<>(); - List sectColors = new ArrayList<>(); - Map subProxyIDs = new HashMap<>(); - for (FaultSection sect : proxyModule.getProxySects()) { - if (sect.getParentSectionName().startsWith(parentName)) { - subProxyIDs.put(sect.getSectionId(), combFaultSects.size()); - combFaultSects.add(sect); + String prefix = parentName.replaceAll("\\W+", "_"); + while (parentName.startsWith("_")) + parentName = parentName.substring(1); + while (parentName.endsWith("_")) + parentName = parentName.substring(0, parentName.length()-1); + while (parentName.contains("__")) + parentName = parentName.replace("__", "_"); + + List proxySects = new ArrayList<>(); + for (FaultSection sect : rupSet.getFaultSectionDataList()) + if (sect.getParentSectionId() == parentID) + proxySects.add(sect); + + int rupID = -1; + double maxRupMag = 0d; + List allRups = new ArrayList<>(rupSet.getRupturesForParentSection(parentID)); + Collections.reverse(allRups); + for (int rupIndex : allRups) { + double mag = rupSet.getMagForRup(rupIndex); + if ((float)mag < (float)maxRupMag) + continue; + if (rupSet.getLengthForRup(rupIndex) > lengthLimit*1e3) // km to m + continue; + List rupSects = rupSet.getFaultSectionDataForRupture(rupIndex); +// if (rupSects.size() != proxySects.size()) +// continue; + boolean allMatch = true; + for (FaultSection sect : rupSects) { + if (sect.getParentSectionId() != parentID) { + allMatch = false; + break; } } - combFaultSects.addAll(proxySects); - mapMaker.setFaultSections(combFaultSects); - for (int i=0; i> rupProxyIndexes = proxyModule.getRupProxySectIndexes(rupID); - for (int i=0; i combFaultSects = new ArrayList<>(); - List sectColors = new ArrayList<>(); - mapMaker.setSkipNaNs(true); - - double len = rup.getAveLength(); - double dip = rup.getAveDip(); - double strike = rup.getAveStrike(); - for (int i=0; i= 0, "No matching rup found for %s. %s; tried %s on parent.", + parentID, parentName, allRups.size()); + + System.out.println(parentName); + System.out.println("Using an M"+(float)maxRupMag+" rupture, len="+(float)rupSet.getLengthForRup(rupID)*1e-3); + + double rupLen = rupSet.getLengthForRup(rupID)*1e-3; + double regBuffer = Math.min(60d, Double.max(30d, rupLen+5d)); + lowerLeft = LocationUtils.location(lowerLeft, 5d*Math.PI/4d, regBuffer); + upperRight = LocationUtils.location(upperRight, Math.PI/4d, regBuffer); + + Region region = new Region(lowerLeft, upperRight); + + UniqueRupture rupUnique = UniqueRupture.forIDs(rupSet.getSectionsIndicesForRup(rupID)); + + LocationList gridLocs = new LocationList(); + for (FaultSection sect : proxySects) { + if (rupUnique.contains(sect.getSectionId())) { + GriddedRegion grid = new GriddedRegion(sect.getZonePolygon(), 0.05d, GriddedRegion.ANCHOR_0_0); + gridLocs.addAll(grid.getNodeList()); } - combFaultSects.addAll(proxySects); - while (sectColors.size() < combFaultSects.size()) - sectColors.add(null); - mapMaker.setSkipNaNs(true); - mapMaker.setFaultSections(combFaultSects); - mapMaker.plotSectColors(sectColors); - mapMaker.plot(outputDiur, "zone_rups_finite", "Finite Fault Ruptures"); - } else if (type == RupType.RAND_STRIKE_FINITE) { - List combFaultSects = new ArrayList<>(); - List sectColors = new ArrayList<>(); - mapMaker.setSkipNaNs(true); - - double len = rup.getAveLength(); - double dip = rup.getAveDip(); - double strike = rup.getAveStrike(); + } + + RuptureSurface rup = rupSet.getSurfaceForRupture(rupID, 1d); + + for (RupType type : types) { + GeographicMapMaker mapMaker = new GeographicMapMaker(region); + mapMaker.setWriteGeoJSON(false); + mapMaker.setScalarThickness(3f); - Random rand = new Random(proxySects.size()); - for (int i=0; i sectColors = new ArrayList<>(); + for (int i=0; i chars = new ArrayList<>(); + for (int i=0; i combFaultSects = new ArrayList<>(); + List sectColors = new ArrayList<>(); + Map subProxyIDs = new HashMap<>(); + for (FaultSection sect : proxyModule.getProxySects()) { + if (sect.getParentSectionName().startsWith(parentName)) { + subProxyIDs.put(sect.getSectionId(), combFaultSects.size()); + combFaultSects.add(sect); + } + } + combFaultSects.addAll(proxySects); + mapMaker.setFaultSections(combFaultSects); + for (int i=0; i> rupProxyIndexes = proxyModule.getRupProxySectIndexes(rupID); + for (int i=0; i combFaultSects = new ArrayList<>(); + List sectColors = new ArrayList<>(); + mapMaker.setSkipNaNs(true); + + double len = rup.getAveLength(); + double dip = rup.getAveDip(); + double strike = rup.getAveStrike(); + for (int i=0; i combFaultSects = new ArrayList<>(); + List sectColors = new ArrayList<>(); + mapMaker.setSkipNaNs(true); + + double len = rup.getAveLength(); + double dip = rup.getAveDip(); + double strike = rup.getAveStrike(); + + Random rand = new Random(proxySects.size()); + for (int i=0; i origSubSects = PRVI25_CrustalDeformationModels.GEOLOGIC.build(fm); + List targetSubSects = PRVI25_CrustalDeformationModels.GEOLOGIC_DIST_AVG.build(fm); + System.out.println("Model has "+numTotSects+" sections ("+(numTotSects-numProxySects)+" regular, " + +numProxySects+" proxy) and "+origSubSects.size()+" subsections"); + GeographicMapMaker mapMaker = new GeographicMapMaker(CRUSTAL_FAULT_MAP_REG, origSubSects); + mapMaker.setWriteGeoJSON(false); + + CPT slipCPT = GMT_CPT_Files.SEQUENTIAL_BATLOW_UNIFORM.instance().rescale(0d, 11d); + slipCPT.setPreferredTickInterval(2d); + + for (boolean target : new boolean[] {false,true}) { + String prefix = target ? "crustal_target_slip_rates" : "crustal_original_slip_rates"; + List subSects = target ? targetSubSects : origSubSects; + List slipRates = getDMSlipRates(subSects, false); + System.out.println("Max crustal is "+slipRates.stream().mapToDouble(D->D).max().getAsDouble()); + + String label = (target ? "Target " : "")+"Slip Rate (mm/yr)"; + mapMaker.plotSectScalars(slipRates, slipCPT, label); + + mapMaker.plot(outputDir, prefix, " "); + mapMaker.setCPTLocation(RectangleEdge.TOP); + PlotSpec linearPlot = mapMaker.buildPlot(" "); + mapMaker.setCPTLocation(RectangleEdge.BOTTOM); + + CPT logSlipCPT = GMT_CPT_Files.SEQUENTIAL_LAJOLLA_UNIFORM.instance().rescale(-1, 1); +// slipCPT.setPreferredTickInterval(2d); + + List logSlipRates = getDMSlipRates(subSects, true); + + mapMaker.plotSectScalars(logSlipRates, logSlipCPT, "Log10 "+label); + + mapMaker.plot(outputDir, prefix+"_log", " "); + PlotSpec logPlot = mapMaker.buildPlot(" "); + + HeadlessGraphPanel gp = PlotUtils.initHeadless(); + + gp.drawGraphPanel(List.of(linearPlot, logPlot), false, false, List.of(mapMaker.getXRange()), + List.of(mapMaker.getYRange(), mapMaker.getYRange())); + + PlotUtils.writePlots(outputDir, prefix+"_combined", gp, mapMaker.getDefaultPlotWidth(), true, true, true, false); + } + + FaultSystemRupSet rupSet = sol.getRupSet(); + List targetSlips = new ArrayList<>(); + Map> targetSlipsByParent = new HashMap<>(); + SectSlipRates rupSetTargets = rupSet.getSectSlipRates(); + for (int i=0; i()); + targetSlipsByParent.get(parentID).add(slip); + } + Map targetAverageSlipsByParent = targetSlipsByParent.entrySet().stream() + .collect(Collectors.toMap(Map.Entry::getKey, + entry -> entry.getValue().stream().mapToDouble(Double::doubleValue).average().orElse(0.0) + )); + + for (boolean includeClassic : new boolean[] {false,true}) { + String prefix = "crustal_solution_slip_rates"; + String solName; + SolutionSlipRates solSlipsModule; + if (includeClassic) { + solSlipsModule = sol.requireModule(SolutionSlipRates.class); + solName = "Solution"; + } else { + solSlipsModule = solNoClassic.requireModule(SolutionSlipRates.class); + prefix += "_no_classic"; + solName = "Solution (excl. classic)"; + } + List solSlips = new ArrayList<>(); + for (int i=0; i> solSlipsByParent = new HashMap<>(); + for (FaultSection sect : rupSet.getFaultSectionDataList()) { + int parentID = sect.getParentSectionId(); + if (!solSlipsByParent.containsKey(parentID)) + solSlipsByParent.put(parentID, new ArrayList<>()); + solSlipsByParent.get(parentID).add(solSlipsModule.get(sect.getSectionId())*1e3); + } + Map solAverageSlipsByParent = solSlipsByParent.entrySet().stream() + .collect(Collectors.toMap(Map.Entry::getKey, + entry -> entry.getValue().stream().mapToDouble(Double::doubleValue).average().orElse(0.0) + )); +// List solSlips = new ArrayList<>(); +// for (FaultSection sect : rupSet.getFaultSectionDataList()) { +// int parentID = sect.getParentSectionId(); +// solSlips.add(solAverageSlipsByParent.get(parentID)); +// if (sect.getSectionId() > 0 && parentID != rupSet.getFaultSectionData(sect.getSectionId()-1).getParentSectionId()) +// System.out.println(sect.getParentSectionName()+": "+solAverageSlipsByParent.get(parentID)); +// } + + mapMaker.plotSectScalars(solSlips, slipCPT, solName+" Slip Rates (mm/yr)"); + mapMaker.plot(outputDir, prefix, " "); + + CPT diffCPT = GMT_CPT_Files.DIVERGING_BAM_UNIFORM.instance().reverse().rescale(-1d, 1d); + CPT pDiffCPT = GMT_CPT_Files.DIVERGING_VIK_UNIFORM.instance().rescale(-20d, 20d); + + mapMaker.plotSectScalars(difference(targetSlips, solSlips), diffCPT, solName+" - Target Slip Rates (mm/yr)"); + mapMaker.plot(outputDir, prefix+"_diff", " "); + + mapMaker.plotSectScalars(percentDifference(targetSlips, solSlips), pDiffCPT, solName+" vs Target Slip Rates (% Difference)"); + mapMaker.plot(outputDir, prefix+"_pDiff", " "); + + plotScatter(outputDir, prefix+"_scatter", targetSlips, solSlips, targetAverageSlipsByParent, solAverageSlipsByParent); + } + } + + private static FaultSystemSolution buildNoClassic(File nodeSolDir) throws IOException { + BranchAverageSolutionCreator baCreator = new BranchAverageSolutionCreator(new BranchWeightProvider.CurrentWeights()); + List> levels = new ArrayList<>(); + levels.add(PRVI25_LogicTreeBranch.SEG); + for (NSHM23_SegmentationModels segModel : NSHM23_SegmentationModels.values()) { + if (segModel.getNodeWeight(null) == 0d || segModel == NSHM23_SegmentationModels.CLASSIC) + continue; + File solFile = new File(nodeSolDir, "SegModel_"+segModel.getFilePrefix()+".zip"); + LogicTreeBranch branch = new LogicTreeBranch<>(levels); + branch.setValue(segModel); + FaultSystemSolution sol = FaultSystemSolution.load(solFile); + baCreator.addSolution(sol, branch); + } + return baCreator.build(); + } + + private static List getDMSlipRates(List subSects, boolean log) { + List slipRates = new ArrayList<>(); + + for (FaultSection sect : subSects) { + double slip = sect.getReducedAveSlipRate(); + if (log) + slip = Math.log10(slip); + slipRates.add(slip); + } + + return slipRates; + } + + private static List difference(List target, List solution) { + List ret = new ArrayList<>(); + for (int i=0; i percentDifference(List target, List solution) { + List ret = new ArrayList<>(); + for (int i=0; i targetSlips, List solSlips, + Map targetParentSlips, Map solParentSlips) throws IOException { + List funcs = new ArrayList<>(); + List chars = new ArrayList<>(); + + Range linearRange = new Range(0d, 13d); + Range logRange = new Range(1e-1, 2e1); + + DefaultXY_DataSet oneToOne = new DefaultXY_DataSet(); + oneToOne.set(linearRange.getLowerBound(), linearRange.getLowerBound()); + oneToOne.set(logRange.getLowerBound(), logRange.getLowerBound()); + oneToOne.set(linearRange.getUpperBound(), linearRange.getUpperBound()); + oneToOne.set(logRange.getUpperBound(), logRange.getUpperBound()); + + funcs.add(oneToOne); + chars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 3f, Color.GRAY)); + + if (solSlips != null) { + DefaultXY_DataSet scatter = new DefaultXY_DataSet(); + for (int i=0; i