From 3a64aecb800162be16b99f4130618bde4242f3e6 Mon Sep 17 00:00:00 2001 From: Kevin Milner Date: Fri, 20 Dec 2024 13:53:37 -0800 Subject: [PATCH] PRVI paper figures and calcs --- ..._LogicTreeInversionRunnerScriptWriter.java | 8 +- .../figures/WUS_HazardChangePageGen.java | 8 +- .../nshm23/wrapper/RewriteExternalHazard.java | 23 +- .../CrustalSubductionLogicTreeCombine.java | 3 +- .../kevin/prvi25/GMMLogicTreeWriter.java | 131 +++--- .../kevin/prvi25/HazardComp2003PageGen.java | 443 ++++++++++++++++++ .../PRVI_SubductionTestRupSetBuilder.java | 3 +- .../kevin/prvi25/SubductionDefModConvert.java | 3 +- ...onInterfaceSlabHazardLogicTreeCombine.java | 7 +- .../prvi25/figures/GridCubeRatePlot.java | 44 +- .../prvi25/figures/IndividualMFDPlots.java | 77 ++- .../figures/PRVI_SubductionSubSectPlots.java | 51 +- .../ProxyFaultRepresentationFigures.java | 2 +- .../prvi25/figures/RupSetStatsTexWriter.java | 10 + .../kevin/prvi25/figures/SlipRateFigures.java | 204 +++++++- .../scratch/kevin/ucerf3/PureScratch.java | 37 +- 16 files changed, 949 insertions(+), 105 deletions(-) create mode 100644 src/main/java/scratch/kevin/prvi25/HazardComp2003PageGen.java diff --git a/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java b/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java index 1977c491..c65e85fd 100644 --- a/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java +++ b/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java @@ -118,6 +118,7 @@ public static void main(String[] args) throws IOException { boolean hazardGridded = false; boolean forceRequiredNonzeroWeight = false; Double forceHazardGridSpacing = null; + long randSeed = 12345678l; File remoteMainDir = new File("/project/scec_608/kmilner/nshm23/batch_inversions"); int remoteTotalThreads = 20; @@ -638,7 +639,11 @@ public static void main(String[] args) throws IOException { levels.remove(i); Preconditions.checkState(levels.size() == origNumLevels -1); individualRandomLevels.add(new PRVI25_CrustalRandomlySampledDeformationModelLevel()); - samplingBranchCountMultiplier = 5; // 5 for each branch +// samplingBranchCountMultiplier = 5; // 5 for each branch +// samplingBranchCountMultiplier = 10; // 10 for each branch +// samplingBranchCountMultiplier = 20; // 20 for each branch + samplingBranchCountMultiplier = 50; // 50 for each branch + randSeed *= samplingBranchCountMultiplier; dirName += "-dmSample"; if (samplingBranchCountMultiplier > 1) dirName += samplingBranchCountMultiplier+"x"; @@ -756,7 +761,6 @@ public static void main(String[] args) throws IOException { // int numSamples = nodes*5; // int numSamples = nodes*4; // long randSeed = System.currentTimeMillis(); - long randSeed = 12345678l; int numSamples = 0; // int numSamples = 450; // int numSamples = 36*10; diff --git a/src/main/java/scratch/kevin/nshm23/figures/WUS_HazardChangePageGen.java b/src/main/java/scratch/kevin/nshm23/figures/WUS_HazardChangePageGen.java index 6de30eb3..40c2dabe 100644 --- a/src/main/java/scratch/kevin/nshm23/figures/WUS_HazardChangePageGen.java +++ b/src/main/java/scratch/kevin/nshm23/figures/WUS_HazardChangePageGen.java @@ -1562,13 +1562,13 @@ else if (hazDiff < 0 && moDiff < 0) private static final DecimalFormat twoDigits = new DecimalFormat("0.00"); private static final DecimalFormat pDF = new DecimalFormat("0.00%"); - static GriddedGeoDataSet asLog10(GriddedGeoDataSet xyz) { + public static GriddedGeoDataSet asLog10(GriddedGeoDataSet xyz) { xyz = xyz.copy(); xyz.log10(); return xyz; } - static GriddedGeoDataSet sumMap(GriddedGeoDataSet map1, GriddedGeoDataSet map2) { + public static GriddedGeoDataSet sumMap(GriddedGeoDataSet map1, GriddedGeoDataSet map2) { Preconditions.checkState(map1.size() == map2.size()); GriddedGeoDataSet ret = new GriddedGeoDataSet(map1.getRegion()); for (int i=0; i tree, Map, LogicTre Map nodeRemaps, String nameAdd, String shortNameAdd) { for (LogicTreeLevel level : tree.getLevels()) { String name = level.getName(); - if (name.toLowerCase().contains("crustal") || name.toLowerCase().contains("subduction")) { + if (name.toLowerCase().contains("crustal") || name.toLowerCase().contains("subduction") + || name.toLowerCase().contains("interface") || name.toLowerCase().contains("slab")) { // keep it as is levelRemaps.put(level, level); for (LogicTreeNode node : level.getNodes()) diff --git a/src/main/java/scratch/kevin/prvi25/GMMLogicTreeWriter.java b/src/main/java/scratch/kevin/prvi25/GMMLogicTreeWriter.java index c4e64518..798c0ef7 100644 --- a/src/main/java/scratch/kevin/prvi25/GMMLogicTreeWriter.java +++ b/src/main/java/scratch/kevin/prvi25/GMMLogicTreeWriter.java @@ -23,6 +23,7 @@ import org.opensha.commons.logicTree.LogicTreeNode; import org.opensha.commons.logicTree.LogicTreeNode.FileBackedNode; import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; +import org.opensha.sha.earthquake.faultSysSolution.hazard.AbstractLogicTreeHazardCombiner; import org.opensha.sha.earthquake.faultSysSolution.hazard.mpj.MPJ_LogicTreeHazardCalc; import org.opensha.sha.earthquake.faultSysSolution.hazard.mpj.MPJ_SiteLogicTreeHazardCurveCalc; import org.opensha.sha.earthquake.faultSysSolution.modules.SolutionLogicTree; @@ -57,6 +58,8 @@ public static void main(String[] args) throws IOException { vs30 = 760d; dirSuffix = "-vs760"; Double sigmaTrunc = 3d; boolean supersample = true; + int erfSamples = -1; + int gmmSamplesPerERF = -1; double[] periods = { 0d, 0.2d, 1d, 5d }; /* @@ -74,30 +77,37 @@ public static void main(String[] args) throws IOException { // // including gridded // int mins = 1440*5; // File sourceTreeFile = new File(sourceDir, "logic_tree_full_gridded.json"); +// erfSamples = 10000; gmmSamplesPerERF = 1; jobSuffix = "_sampled"; logicTreeOutputName = "logic_tree_full_gridded_sampled.json"; //// File sourceTreeFile = new File(sourceDir, "logic_tree_full_gridded_sampled.json"); jobSuffix = "_sampled"; // IncludeBackgroundOption bgOp = IncludeBackgroundOption.INCLUDE; /* - * Interface + * Interface separate slab and interface * - * do supra-seis, then gridded-only, then combine + * do supra-seis, then each gridded-only, then interface combine + * + * then need to separately combine the slab logic tree */ -// List> gmmLevels = PRVI25_LogicTreeBranch.levelsInterfaceGMM; -// File sourceDir = SUBDUCTION_DIR; -// File outputDir = new File(sourceDir.getParentFile(), sourceDir.getName()+"-gmTreeCalcs"+dirSuffix); -// // supra-seis only -//// File sourceTreeFile = new File(sourceDir, "logic_tree.json"); -//// int mins = 1440; -//// IncludeBackgroundOption bgOp = IncludeBackgroundOption.EXCLUDE; -// // interface gridded only -//// int mins = 1440; -//// File sourceTreeFile = new File(sourceDir, "logic_tree_gridded_only.json"); -//// logicTreeOutputName = "logic_tree_gridded_interface_only.json"; -//// IncludeBackgroundOption bgOp = IncludeBackgroundOption.ONLY; -//// forceInputFileName = "results_gridded_branches_interface_only.zip"; -//// jobSuffix = "_interface"; -//// outputSuffix = jobSuffix; -// // interface both (combine only) + List> gmmLevels = PRVI25_LogicTreeBranch.levelsInterfaceGMM; + File sourceDir = SUBDUCTION_DIR; + File outputDir = new File(sourceDir.getParentFile(), sourceDir.getName()+"-gmTreeCalcs"+dirSuffix); + // supra-seis only +// File sourceTreeFile = new File(sourceDir, "logic_tree.json"); +// int mins = 1440; +// IncludeBackgroundOption bgOp = IncludeBackgroundOption.EXCLUDE; + // interface gridded only + int mins = 1440; +// File sourceTreeFile = new File(sourceDir, "logic_tree_gridded_only.json"); + File sourceTreeFile = new File(sourceDir, "logic_tree_full_gridded_for_only_calc.json"); + logicTreeOutputName = "logic_tree_gridded_interface_only.json"; + IncludeBackgroundOption bgOp = IncludeBackgroundOption.ONLY; + // this was for if gridded only depended on FM but it also depends on scale +// forceInputFileName = "results_gridded_branches_interface_only.zip"; + // use this one because it has scaling relationship specific gridded models + forceInputFileName = "results_full_gridded_interface_only.zip"; + jobSuffix = "_interface"; + outputSuffix = jobSuffix; + // interface both (combine only) // combineOnly = true; // int mins = 1440; // forceInputFileName = "results_full_gridded_interface_only.zip"; @@ -125,26 +135,26 @@ public static void main(String[] args) throws IOException { /* * Branch averaged (GMM-only) */ - List> gmmLevels = PRVI25_LogicTreeBranch.levelsCombinedGMM; - File sourceDir = COMBINED_DIR; - File outputDir = new File(sourceDir.getParentFile(), sourceDir.getName()+"-ba_only-gmTreeCalcs"+dirSuffix); - // write out a SLT that only contains that node - File sourceTreeFile = new File(outputDir, "fake_erf_logic_tree.json"); - FileBackedLevel fakeLevel = new FileBackedLevel("ERF Model", "ERF", - List.of(new FileBackedNode("Branch Averaged ERF", "BranchAveragedERF", 1d, "BA_ERF"))); - LogicTree tempTree = LogicTree.buildExhaustive(List.of(fakeLevel), true); - Preconditions.checkState(tempTree.size() == 1); - File sourceFile = new File(outputDir, "fake_erf_slt.zip"); - SolutionLogicTree.FileBuilder builder = new SolutionLogicTree.FileBuilder(sourceFile); - builder.setSerializeGridded(true); - builder.solution(FaultSystemSolution.load(COMBINED_SOL), tempTree.getBranch(0)); - builder.close(); - forceInputFileName = sourceFile.getName(); - tempTree.write(sourceTreeFile); - logicTreeOutputName = "logic_tree.json"; - sourceDir = outputDir; - int mins = 1440; - IncludeBackgroundOption bgOp = IncludeBackgroundOption.INCLUDE; +// List> gmmLevels = PRVI25_LogicTreeBranch.levelsCombinedGMM; +// File sourceDir = COMBINED_DIR; +// File outputDir = new File(sourceDir.getParentFile(), sourceDir.getName()+"-ba_only-gmTreeCalcs"+dirSuffix); +// // write out a SLT that only contains that node +// File sourceTreeFile = new File(outputDir, "fake_erf_logic_tree.json"); +// FileBackedLevel fakeLevel = new FileBackedLevel("ERF Model", "ERF", +// List.of(new FileBackedNode("Branch Averaged ERF", "BranchAveragedERF", 1d, "BA_ERF"))); +// LogicTree tempTree = LogicTree.buildExhaustive(List.of(fakeLevel), true); +// Preconditions.checkState(tempTree.size() == 1); +// File sourceFile = new File(outputDir, "fake_erf_slt.zip"); +// SolutionLogicTree.FileBuilder builder = new SolutionLogicTree.FileBuilder(sourceFile); +// builder.setSerializeGridded(true); +// builder.solution(FaultSystemSolution.load(COMBINED_SOL), tempTree.getBranch(0)); +// builder.close(); +// forceInputFileName = sourceFile.getName(); +// tempTree.write(sourceTreeFile); +// logicTreeOutputName = "logic_tree.json"; +// sourceDir = outputDir; +// int mins = 1440; +// IncludeBackgroundOption bgOp = IncludeBackgroundOption.INCLUDE; // FOR ALL @@ -159,23 +169,32 @@ public static void main(String[] args) throws IOException { File gridRegFile = new File(outputDir, "gridded_region.geojson"); Feature.write(gridReg.toFeature(), gridRegFile); - List> combLevels = new ArrayList<>(); - combLevels.addAll(erfTree.getLevels()); - combLevels.addAll(gmmLevels); - - List> combBranches = new ArrayList<>(erfTree.size()*gmmTree.size()); - - for (LogicTreeBranch branch : erfTree) { - for (LogicTreeBranch gmmBranch : gmmTree) { - LogicTreeBranch combBranch = new LogicTreeBranch<>(combLevels); - for (LogicTreeNode node : branch) - combBranch.setValue(node); - for (LogicTreeNode node : gmmBranch) - combBranch.setValue(node); - combBranches.add(combBranch); + LogicTree logicTree; + if (gmmSamplesPerERF > 0) { + System.out.println("Pairwise sampling with pairwise="+gmmSamplesPerERF+" branches"); + if (erfSamples <= 0) + erfSamples = erfTree.size(); + logicTree = AbstractLogicTreeHazardCombiner.pairwiseSampleLogicTrees(erfTree, gmmTree, erfSamples, gmmSamplesPerERF); + } else { + List> combLevels = new ArrayList<>(); + combLevels.addAll(erfTree.getLevels()); + combLevels.addAll(gmmLevels); + + List> combBranches = new ArrayList<>(erfTree.size()*gmmTree.size()); + + for (LogicTreeBranch branch : erfTree) { + for (LogicTreeBranch gmmBranch : gmmTree) { + LogicTreeBranch combBranch = new LogicTreeBranch<>(combLevels); + for (LogicTreeNode node : branch) + combBranch.setValue(node); + for (LogicTreeNode node : gmmBranch) + combBranch.setValue(node); + combBranches.add(combBranch); + } } + + logicTree = LogicTree.fromExisting(combLevels, combBranches); } - LogicTree logicTree = LogicTree.fromExisting(combLevels, combBranches); File localLogicTree = new File(outputDir, logicTreeOutputName == null ? sourceTreeFile.getName() : logicTreeOutputName); logicTree.write(localLogicTree); @@ -246,7 +265,7 @@ else if (bgOp == IncludeBackgroundOption.EXCLUDE) argz += " --gridded-seis "+bgOp.name(); String logicTreePath = dirPath+"/"+localLogicTree.getName(); argz += " --logic-tree "+logicTreePath; - if (bgOp == IncludeBackgroundOption.ONLY || bgOp == IncludeBackgroundOption.INCLUDE) + if (!inputFileName.contains("interface_only") && (bgOp == IncludeBackgroundOption.ONLY || bgOp == IncludeBackgroundOption.INCLUDE)) argz += " --quick-grid-calc"; if (combineOnly) argz += " --combine-only"; @@ -254,7 +273,7 @@ else if (bgOp == IncludeBackgroundOption.EXCLUDE) if (vs30 != null) argz += " --vs30 "+vs30.floatValue(); if (supersample) - argz += " --supersample"; + argz += " --supersample-quick"; if (sigmaTrunc != null) argz += " --gmm-sigma-trunc-one-sided "+sigmaTrunc.floatValue(); if (periods != null) { @@ -292,7 +311,7 @@ else if (bgOp == IncludeBackgroundOption.EXCLUDE) if (vs30 != null) argz += " --vs30 "+vs30.floatValue(); if (supersample) - argz += " --supersample"; + argz += " --supersample-quick"; if (sigmaTrunc != null) argz += " --gmm-sigma-trunc-one-sided "+sigmaTrunc.floatValue(); if (periods != null) { diff --git a/src/main/java/scratch/kevin/prvi25/HazardComp2003PageGen.java b/src/main/java/scratch/kevin/prvi25/HazardComp2003PageGen.java new file mode 100644 index 00000000..bb1969ef --- /dev/null +++ b/src/main/java/scratch/kevin/prvi25/HazardComp2003PageGen.java @@ -0,0 +1,443 @@ +package scratch.kevin.prvi25; + +import java.awt.Color; +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; +import java.util.List; + +import org.opensha.commons.data.CSVFile; +import org.opensha.commons.data.function.DiscretizedFunc; +import org.opensha.commons.data.function.LightFixedXFunc; +import org.opensha.commons.data.xyz.GriddedGeoDataSet; +import org.opensha.commons.geo.GriddedRegion; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.LocationUtils; +import org.opensha.commons.gui.plot.GeographicMapMaker; +import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; +import org.opensha.commons.util.MarkdownUtils; +import org.opensha.commons.util.MarkdownUtils.TableBuilder; +import org.opensha.commons.util.cpt.CPT; +import org.opensha.sha.earthquake.faultSysSolution.hazard.mpj.MPJ_LogicTreeHazardCalc; +import org.opensha.sha.earthquake.faultSysSolution.ruptures.util.RupSetMapMaker; +import org.opensha.sha.earthquake.faultSysSolution.util.SolHazardMapCalc; +import org.opensha.sha.earthquake.faultSysSolution.util.SolHazardMapCalc.ReturnPeriods; +import org.opensha.sha.earthquake.param.IncludeBackgroundOption; +import org.opensha.sha.earthquake.rupForecastImpl.nshm23.util.NSHM23_RegionLoader; +import org.opensha.sha.earthquake.rupForecastImpl.prvi25.util.PRVI25_RegionLoader; + +import com.google.common.base.Preconditions; + +import scratch.kevin.nshm23.figures.MethodsAndIngredientsHazChangeFigures; +import scratch.kevin.nshm23.figures.WUS_HazardChangePageGen; +import scratch.kevin.prvi25.figures.PRVI_Paths; + +public class HazardComp2003PageGen { + + private static enum CompType { + TOTAL, + CRUSTAL_FAULTS, + CRUSTAL_GRIDDED, + CRUSTAL, + INTERFACE, + SLAB + }; + + public static void main(String[] args) throws IOException { + String imtName = "PGA"; + String imtDir = "PGA"; + double period = 0d; + +// String imtName = "1s SA"; +// String imtDir = "SA1P0"; +// double period = 1d; + +// ReturnPeriods[] rps = ReturnPeriods.values(); + ReturnPeriods[] rps = { ReturnPeriods.TWO_IN_50, ReturnPeriods.TEN_IN_50 }; + + String name25 = "NSHM25"; + String name03 = "NSHM03"; + + File inputDirFull25 = new File(PRVI_Paths.INV_DIR, PRVI_Paths.COMBINED_DIR.getName()+"-ba_only-vs760"); + File inputDirCrustal25 = new File(PRVI_Paths.INV_DIR, PRVI_Paths.CRUSTAL_DIR.getName()+"-ba_only-vs760"); + File inputDirInterface25 = new File(PRVI_Paths.INV_DIR, PRVI_Paths.SUBDUCTION_DIR.getName()+"-ba_only-INTERFACE_only-vs760"); + File inputDirSlab25 = new File(PRVI_Paths.INV_DIR, PRVI_Paths.SUBDUCTION_DIR.getName()+"-ba_only-SLAB_only-vs760"); + File outputDir = new File(PRVI_Paths.COMBINED_DIR, "nshm03_erf_comparisons_"+imtDir); + double sourceSpacing = 0.01; + File sourcesDir03 = new File("/home/kevin/OpenSHA/nshm23/nshmp-haz-models/ext_hazard_calcs/" + + "prvi-t2-2003-ERF-2025-vB1-GMMs-2025conf-0p01-vs760-20241213-49f58ecb02d600/vs30-760/"+imtDir+"/source"); + + System.out.println("Output dir: "+outputDir.getAbsolutePath()); + Preconditions.checkState(outputDir.exists() || outputDir.mkdir()); + File resourcesDir = new File(outputDir, "resources"); + Preconditions.checkState(resourcesDir.exists() || resourcesDir.mkdir()); + + GriddedRegion mapReg = new GriddedRegion(PRVI25_RegionLoader.loadPRVI_MapExtents(), + sourceSpacing, GriddedRegion.ANCHOR_0_0); + + DiscretizedFunc[] crustalFault25 = loadRegularCurves(inputDirCrustal25, IncludeBackgroundOption.EXCLUDE, mapReg, period); + DiscretizedFunc[] crustalGrid25 = loadRegularCurves(inputDirCrustal25, IncludeBackgroundOption.ONLY, mapReg, period); + DiscretizedFunc[] interface25 = loadRegularCurves(inputDirInterface25, IncludeBackgroundOption.INCLUDE, mapReg, period); + DiscretizedFunc[] slab25 = loadRegularCurves(inputDirSlab25, IncludeBackgroundOption.ONLY, mapReg, period); + DiscretizedFunc[] full25 = loadRegularCurves(inputDirFull25, IncludeBackgroundOption.INCLUDE, mapReg, period); + + DiscretizedFunc[] crustalFault03 = add( + loadExtCurves(new File(sourcesDir03, "FAULT/curves.csv"), mapReg), + loadExtCurves(new File(sourcesDir03, "ZONE/curves.csv"), mapReg) + ); + DiscretizedFunc[] crustalGrid03 = loadExtCurves(new File(sourcesDir03, "GRID/curves.csv"), mapReg); + DiscretizedFunc[] interface03 = loadExtCurves(new File(sourcesDir03, "INTERFACE/curves.csv"), mapReg); + DiscretizedFunc[] slab03 = loadExtCurves(new File(sourcesDir03, "SLAB/curves.csv"), mapReg); + DiscretizedFunc[] full03 = add(crustalFault03, crustalGrid03, interface03, slab03); + + // convert to probabilities + crustalFault03 = ratesToProbs(crustalFault03); + crustalGrid03 = ratesToProbs(crustalGrid03); + interface03 = ratesToProbs(interface03); + slab03 = ratesToProbs(slab03); + full03 = ratesToProbs(full03); + + // interpolate 03 onto our x-value spacing + crustalFault03 = interpolateToMatch(crustalFault03, full25[0]); + crustalGrid03 = interpolateToMatch(crustalGrid03, full25[0]); + interface03 = interpolateToMatch(interface03, full25[0]); + slab03 = interpolateToMatch(slab03, full25[0]); + full03 = interpolateToMatch(full03, full25[0]); + + // these maps substitute a single nshm03 ingredient at a time, keeping everything else at 25 + DiscretizedFunc[] withCrustalFaults03 = add(crustalFault03, crustalGrid25, interface25, slab25); + DiscretizedFunc[] withCrustalGrid03 = add(crustalFault25, crustalGrid03, interface25, slab25); + DiscretizedFunc[] withCrustal03 = add(crustalFault03, crustalGrid03, interface25, slab25); + DiscretizedFunc[] withInterface03 = add(crustalFault25, crustalGrid25, interface03, slab25); + DiscretizedFunc[] withSlab03 = add(crustalFault25, crustalGrid25, interface25, slab03); + + GeographicMapMaker mapMaker = new RupSetMapMaker(List.of(), mapReg); +// mapMaker.setDefaultPlotWidth(1000); + + Color transparent = new Color(255, 255, 255, 0); + + CPT hazCPT = GMT_CPT_Files.RAINBOW_UNIFORM.instance().rescale(-3, 1); + hazCPT.setNanColor(transparent); + + CPT pDiffCPT = MethodsAndIngredientsHazChangeFigures.getCenterMaskedCPT(GMT_CPT_Files.DIVERGING_VIK_UNIFORM.instance(), 10d, 50d); + pDiffCPT.setNanColor(transparent); + + CPT diffCPT = GMT_CPT_Files.DIVERGING_BAM_UNIFORM.instance().reverse().rescale(-0.2d, 0.2d); + diffCPT.setNanColor(transparent); + + List lines = new ArrayList<>(); + + lines.add("# Hazard Comparisons, "+name25+" vs "+name03); + lines.add(""); + + lines.add("This page compares "+name25+" and "+name03+" hazard maps. " + + "Each individual ERF contributor to hazard changes are plotted separately, but GMMs are held constant."); + lines.add(""); + + String csvStr = "Download map CSVs:"; + for (ReturnPeriods rp : rps) { + System.out.println("Building CSV for "+rp); + CSVFile csv = new CSVFile<>(true); + csv.addLine("Location Index", "Latitude", "Longitude", name25, name03, + name03+" Crustal Faults", name03+" Crustal Gridded", name03+" Crustal", + name03+" Interface", name03+" Slab"); + + GriddedGeoDataSet[] maps = { + curvestoMap(full25, mapReg, rp), + curvestoMap(full03, mapReg, rp), + curvestoMap(withCrustalFaults03, mapReg, rp), + curvestoMap(withCrustalGrid03, mapReg, rp), + curvestoMap(withCrustal03, mapReg, rp), + curvestoMap(withInterface03, mapReg, rp), + curvestoMap(withSlab03, mapReg, rp) + }; + + for (int i=0; i line = new ArrayList<>(csv.getNumCols()); + + line.add(i+""); + Location loc = maps[0].getLocation(i); + line.add((float)loc.getLatitude()+""); + line.add((float)loc.getLongitude()+""); + for (GriddedGeoDataSet map : maps) + line.add((float)map.get(i)+""); + + csv.addLine(line); + } + + File csvFile = new File(resourcesDir, "maps_"+rp.name()+".csv"); + csv.writeToFile(csvFile); + + csvStr += " [__"+csvFile.getName()+"__]("+resourcesDir.getName()+"/"+csvFile.getName()+")"; + } + lines.add(csvStr); + lines.add(""); + + int tocIndex = lines.size(); + String topLink = "*[(top)](#table-of-contents)*"; + + for (CompType type : CompType.values()) { + DiscretizedFunc[] curves25, curves03; + String label, description; + String mapLabelAdd; + + switch (type) { + case TOTAL: + curves25 = full25; + curves03 = full03; + label = "Full Model Comparison"; + description = "Full model comparison, including fault, gridded, and subduction model changes."; + mapLabelAdd = ""; + break; + case CRUSTAL_FAULTS: + curves25 = full25; + curves03 = withCrustalFaults03; + label = "Crustal Fault-Only Comparison"; + description = "Crustal fault-only comparison, holding gridded and subduction sources constant (using those from "+name25+")."; + mapLabelAdd = "Crustal Faults, "; + break; + case CRUSTAL_GRIDDED: + curves25 = full25; + curves03 = withCrustalGrid03; + label = "Crustal Gridded-Only Comparison"; + description = "Crustal Gridded-only comparison, holding crustal fault and subduction sources constant (using those from "+name25+")."; + mapLabelAdd = "Crustal Gridded, "; + break; + case CRUSTAL: + curves25 = full25; + curves03 = withCrustal03; + label = "Crustal-Only (Fault & Gridded) Comparison"; + description = "Crustal-only comparison, holding subduction sources constant (using those from "+name25+")."; + mapLabelAdd = "Crustal, "; + break; + case INTERFACE: + curves25 = full25; + curves03 = withInterface03; + label = "Interface-Only Comparison"; + description = "Subduction interface-only comparison, crustal and slab sources constant (using those from "+name25+")."; + mapLabelAdd = "Interface, "; + break; + case SLAB: + curves25 = full25; + curves03 = withSlab03; + label = "Slab-Only Comparison"; + description = "Subduction intraslab-only comparison, crustal and interface sources constant (using those from "+name25+")."; + mapLabelAdd = "Slab, "; + break; + + default: + throw new IllegalStateException(); + } + + lines.add("## "+label); + lines.add(topLink); lines.add(""); + + lines.add(description); + lines.add(""); + + for (ReturnPeriods rp : rps) { + System.out.println("Plotting "+label+", "+rp.label); + + lines.add("### "+label+", "+rp.label); + lines.add(topLink); lines.add(""); + + String hazLabel = imtName+", "+rp.label; + String prefix = type.name()+"_"+rp.name(); + + GriddedGeoDataSet map25 = curvestoMap(curves25, mapReg, rp); + GriddedGeoDataSet map03 = curvestoMap(curves03, mapReg, rp); + + TableBuilder table = MarkdownUtils.tableBuilder(); + + table.addLine(name25, name03); + + table.initNewLine(); + + mapMaker.plotXYZData(asLog10(map25), hazCPT, mapLabelAdd+name25+", "+hazLabel+" (g)"); + mapMaker.plot(resourcesDir, prefix+"_nshm25", " "); + table.addColumn("![Map]("+resourcesDir.getName()+"/"+prefix+"_nshm25.png)"); + mapMaker.plotXYZData(asLog10(map03), hazCPT, mapLabelAdd+name03+", "+hazLabel+" (g)"); + mapMaker.plot(resourcesDir, prefix+"_nshm03", " "); + table.addColumn("![Map]("+resourcesDir.getName()+"/"+prefix+"_nshm03.png)"); + + table.finalizeLine(); + + table.addLine(MarkdownUtils.boldCentered("Ratio"), MarkdownUtils.boldCentered("Difference")); + + GriddedGeoDataSet pDiff = WUS_HazardChangePageGen.mapPDiff(map25, map03); + GriddedGeoDataSet diff = WUS_HazardChangePageGen.mapDiff(map25, map03); + + table.initNewLine(); + + mapMaker.plotXYZData(pDiff, pDiffCPT, mapLabelAdd+name25+" vs "+name03+", % Change, "+hazLabel); + mapMaker.plot(resourcesDir, prefix+"_pDiff", " "); + table.addColumn("![Map]("+resourcesDir.getName()+"/"+prefix+"_pDiff.png)"); + mapMaker.plotXYZData(diff, diffCPT, mapLabelAdd+name25+" - "+name03+", "+hazLabel+" (g)"); + mapMaker.plot(resourcesDir, prefix+"_diff", " "); + table.addColumn("![Map]("+resourcesDir.getName()+"/"+prefix+"_diff.png)"); + + table.finalizeLine(); + + lines.addAll(table.build()); + lines.add(""); + } + } + + // add TOC + lines.addAll(tocIndex, MarkdownUtils.buildTOC(lines, 2)); + lines.add(tocIndex, "## Table Of Contents"); + + // write markdown + MarkdownUtils.writeReadmeAndHTML(lines, outputDir); + } + + private static DiscretizedFunc[] loadRegularCurves(File runDir, IncludeBackgroundOption bgOp, + GriddedRegion gridReg, double period) throws IOException { + String hazardPrefix = "hazard_"+(float)gridReg.getSpacing()+"deg"; + hazardPrefix += "_grid_seis_"; + hazardPrefix += bgOp.name(); + File resultsDir = new File(runDir, "results"); + File subDir = new File(resultsDir, hazardPrefix); + Preconditions.checkState(subDir.exists(), "Doesn't exist: %s", subDir.getAbsolutePath()); + + File curvesFile = new File(subDir, SolHazardMapCalc.getCSV_FileName("curves", period)); + if (!curvesFile.exists()) + curvesFile = new File(curvesFile.getAbsolutePath()+".gz"); + Preconditions.checkState(curvesFile.exists(), "Doesn't exist: %s", curvesFile.getAbsolutePath()); + + return SolHazardMapCalc.loadCurvesCSV(CSVFile.readFile(curvesFile, true), gridReg); + } + + private static DiscretizedFunc[] loadExtCurves(File csvFile, GriddedRegion gridReg) throws IOException { + System.out.println("Loading external curves from "+csvFile.getAbsolutePath()); + CSVFile csv = CSVFile.readFile(csvFile, true); + + double[] xVals = new double[csv.getNumCols()-2]; + for (int i=0; i curve.getMaxX()) + interpY[i] = 0d; + else + interpY[i] = curve.getInterpolatedY_inLogXDomain(x); + } + ret[n] = new LightFixedXFunc(xVals, interpY); + } + return ret; + } + + private static DiscretizedFunc[] add(DiscretizedFunc[]... allCurves) { + DiscretizedFunc[] ret = new DiscretizedFunc[allCurves[0].length]; + + double[] xVals = new double[allCurves[0][0].size()]; + for (int i=0; i curve.getMaxY()) + val = 0d; + else if (rp.oneYearProb < curve.getMinY()) + // saturated + val = curve.getMaxX(); + else + val = curve.getFirstInterpolatedX_inLogXLogYDomain(rp.oneYearProb); + ret.set(i, val); + } + + return ret; + } + + private static GriddedGeoDataSet asLog10(GriddedGeoDataSet xyz) { + xyz = zerosToNaNs(xyz); + xyz.log10(); + return xyz; + } + + private static GriddedGeoDataSet zerosToNaNs(GriddedGeoDataSet xyz) { + xyz = xyz.copy(); + for (int i=0; i debugPrefixes = HashBasedTable.create(); diff --git a/src/main/java/scratch/kevin/prvi25/SubductionInterfaceSlabHazardLogicTreeCombine.java b/src/main/java/scratch/kevin/prvi25/SubductionInterfaceSlabHazardLogicTreeCombine.java index 2f40b693..c73f3c77 100644 --- a/src/main/java/scratch/kevin/prvi25/SubductionInterfaceSlabHazardLogicTreeCombine.java +++ b/src/main/java/scratch/kevin/prvi25/SubductionInterfaceSlabHazardLogicTreeCombine.java @@ -137,9 +137,10 @@ private static GriddedRegion loadGridReg(File regFile) throws IOException { private static void remapTree(LogicTree tree, Map, LogicTreeLevel> levelRemaps, Map nodeRemaps, String nameAdd, String shortNameAdd) { for (LogicTreeLevel level : tree.getLevels()) { - if (ScalarIMR_ParamsLogicTreeNode.class.isAssignableFrom(level.getType())) { + String name = level.getName(); + if (ScalarIMR_ParamsLogicTreeNode.class.isAssignableFrom(level.getType()) + && !name.toLowerCase().contains("interface") && !name.toLowerCase().contains("slab")) { // need to remap - String name = level.getName(); List modNodes = new ArrayList<>(); for (LogicTreeNode node : level.getNodes()) { FileBackedNode modNode = new FileBackedNode(nameAdd+" "+node.getName(), node.getShortName(), @@ -168,7 +169,7 @@ protected void remapOuterTree(LogicTree tree, Map, LogicTre @Override protected void remapInnerTree(LogicTree tree, Map, LogicTreeLevel> levelRemaps, Map nodeRemaps) { - remapTree(tree, levelRemaps, nodeRemaps, "Slab", ""); + remapTree(tree, levelRemaps, nodeRemaps, "Intraslab", ""); } @Override diff --git a/src/main/java/scratch/kevin/prvi25/figures/GridCubeRatePlot.java b/src/main/java/scratch/kevin/prvi25/figures/GridCubeRatePlot.java index d0f26b40..cf07a183 100644 --- a/src/main/java/scratch/kevin/prvi25/figures/GridCubeRatePlot.java +++ b/src/main/java/scratch/kevin/prvi25/figures/GridCubeRatePlot.java @@ -115,10 +115,10 @@ public static void main(String[] args) throws IOException { props.set(GeoJSONFaultSection.LOW_DEPTH, lowDepth); props.set(GeoJSONFaultSection.UPPER_DEPTH, upDepth); props.set(GeoJSONFaultSection.SLIP_RATE, 10d); - LineString geom = new LineString(traceStart, traceEnd); + Geometry 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))); +// Polygon proxyZone = new Polygon(new Region(new Location(traceStart.lat, -0.5*gridSpacing), +// new Location(traceEnd.lat, 0.5*gridSpacing))); // geom = new GeometryCollection(geom, proxyZone); Feature sectFeature = new Feature(geom, props); GeoJSONFaultSection sect = GeoJSONFaultSection.fromFeature(sectFeature); @@ -256,6 +256,21 @@ public static void main(String[] args) throws IOException { sectXY.set(sectTopLoc.lon*lonSpacingKM, sectTopLoc.depth); sectXY.set(sectBotLoc.lon*lonSpacingKM, sectBotLoc.depth); + XY_DataSet polyXY = null; + if (sect.isProxyFault() && sect.getZonePolygon() != null) { + Region proxyReg = sect.getZonePolygon(); + double left = proxyReg.getMinLon()*lonSpacingKM; + double right = proxyReg.getMaxLon()*lonSpacingKM; + polyXY = new DefaultXY_DataSet(); + polyXY.set(left, sectTopLoc.depth); + polyXY.set(right, sectTopLoc.depth); + polyXY.set(right, sectBotLoc.depth); + polyXY.set(left, sectBotLoc.depth); + polyXY.set(polyXY.get(0)); + + sectXY = null; + } + CPT assocCPT = GMT_CPT_Files.SEQUENTIAL_LAJOLLA_UNIFORM.instance().rescale(0d, 1d); assocCPT.setPreferredTickInterval(0.1); calcXYZValues(assocXYZ, xyzMappedCubes, new Function() { @@ -277,7 +292,7 @@ public Double apply(Integer cubeIndex) { DecimalFormat fractDF = new DecimalFormat("0.00"); XYZPlotSpec assocPlot = buildPlot(assocXYZ, cellXRanges, cellAssocs, new FormatFunc(fractDF), - sectXY, assocCPT, "Fractional Association", " "); + sectXY, polyXY, assocCPT, "Fractional Association", " "); HeadlessGraphPanel gp = PlotUtils.initHeadless(); @@ -337,14 +352,14 @@ public String apply(Double t) { }; XYZPlotSpec cellRatePlot = buildPlot(totalNuclXYZ, cellXRanges, totCellRates, rateFormat, - sectXY, nuclRateCPT, "Gridded M>5 Nucleation Rate", " "); + sectXY, polyXY, nuclRateCPT, "Log10 Gridded M>5 Nucleation Rate", " "); gp.drawGraphPanel(cellRatePlot, false, false, xRange, yRange); PlotUtils.writePlots(outputDir, "cube_nucl_rate_tot", gp, width, false, true, true, false); XYZPlotSpec cellRateM7Plot = buildPlot(m7NuclXYZ, cellXRanges, m7CellRates, rateFormat, - sectXY, nuclRateCPT, "Gridded M>7 Nucleation Rate", " "); + sectXY, polyXY, nuclRateCPT, "Log10 Gridded M>7 Nucleation Rate", " "); gp.drawGraphPanel(cellRateM7Plot, false, false, xRange, yRange); @@ -417,7 +432,7 @@ public String apply(Double t) { } private static XYZPlotSpec buildPlot(EvenlyDiscrXYZ_DataSet xyz, List cellXRanges, List cellValues, - Function cellDF, XY_DataSet sectXY, CPT cpt, String zLabel, String title) { + Function cellDF, XY_DataSet sectXY, XY_DataSet polyXY, CPT cpt, String zLabel, String title) { List funcs = new ArrayList<>(); List chars = new ArrayList<>(); @@ -456,11 +471,18 @@ else if (i == cellValues.size()-1) } } - funcs.add(sectXY); - chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 7f, Color.WHITE)); + if (sectXY != null) { + funcs.add(sectXY); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 7f, Color.WHITE)); + + funcs.add(sectXY); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 4f, Color.BLACK)); + } - funcs.add(sectXY); - chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 4f, Color.BLACK)); + if (polyXY != null) { + funcs.add(polyXY); + chars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 4f, new Color(0, 0, 0, 160))); + } XYZPlotSpec plot = new XYZPlotSpec(xyz, funcs, chars, cpt, title, "X (km)", "Depth (km)", zLabel); plot.setYAxisInverted(true); diff --git a/src/main/java/scratch/kevin/prvi25/figures/IndividualMFDPlots.java b/src/main/java/scratch/kevin/prvi25/figures/IndividualMFDPlots.java index 47632338..ada7389f 100644 --- a/src/main/java/scratch/kevin/prvi25/figures/IndividualMFDPlots.java +++ b/src/main/java/scratch/kevin/prvi25/figures/IndividualMFDPlots.java @@ -163,7 +163,14 @@ else if (trt == TectonicRegionType.SUBDUCTION_INTERFACE) griddedDists = new UncertainBoundedIncrMagFreqDist[] { PRVI25_RegionalSeismicity.getBounded(seisReg, refMFD, PRVI25_GridSourceBuilder.SLAB_MMAX) }; } - if (trt != TectonicRegionType.SUBDUCTION_SLAB) { + if (trt == TectonicRegionType.SUBDUCTION_SLAB) { + if (seisReg == PRVI25_SeismicityRegions.CAR_INTRASLAB) { + // combined slab plot + plotMultiSlab(subductionOutputDir, "subduction_mfds_slab_combined", new Range(5d, 8d), + PRVI25_RegionalSeismicity.getBounded(PRVI25_SeismicityRegions.CAR_INTRASLAB, refMFD, PRVI25_GridSourceBuilder.SLAB_MMAX), + PRVI25_RegionalSeismicity.getBounded(PRVI25_SeismicityRegions.MUE_INTRASLAB, refMFD, PRVI25_GridSourceBuilder.SLAB_MMAX)); + } + } else { if (r > 0) texFW.write("% "+seisReg.name()+"subset "+r+"\n"); else @@ -177,15 +184,34 @@ else if (trt == TectonicRegionType.SUBDUCTION_INTERFACE) texFW.write(LaTeXUtils.defineValueCommand(texPrefix+"ObsMFivePercent", LaTeXUtils.numberAsPercent(100d*obsM5/origM5, 0), false)+"\n"); } + double onFaultTotalRate = onFaultMean.calcSumOfY_Vals(); texFW.write(LaTeXUtils.defineValueCommand(texPrefix+"SupraRate", - LaTeXUtils.numberExpFormatSigFigs(onFaultMean.calcSumOfY_Vals(), 3), false)+"\n"); - double onFaultRI = 1d/onFaultMean.calcSumOfY_Vals(); + LaTeXUtils.numberExpFormatSigFigs(onFaultTotalRate, 3), false)+"\n"); + double onFaultRI = 1d/onFaultTotalRate; if (onFaultRI > 5d) texFW.write(LaTeXUtils.defineValueCommand(texPrefix+"SupraRI", LaTeXUtils.groupedIntNumber(onFaultRI), false)+"\n"); else texFW.write(LaTeXUtils.defineValueCommand(texPrefix+"SupraRI", LaTeXUtils.numberExpFormatFixedDecimal(onFaultRI, 1), false)+"\n"); + if (seisReg == PRVI25_SeismicityRegions.CRUSTAL && r == 0) { + double multiFaultRate = 0d; + for (int rupIndex=0; rupIndex rupSects = sol.getRupSet().getFaultSectionDataForRupture(rupIndex); + for (int i=1; !multiFault&&i funcs = new ArrayList<>(); + List chars = new ArrayList<>(); + + mueDist.setName("Muertos Preferred"); + funcs.add(mueDist); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, mueColor)); + + mueDist.getUpper().setName("Muertos High & Low"); + funcs.add(mueDist.getUpper()); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 1f, mueColor)); + + mueDist.getLower().setName(null); + funcs.add(mueDist.getLower()); + chars.add(chars.get(chars.size()-1)); + + carDist.setName("Caribbean Preferred"); + funcs.add(carDist); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 3f, carColor)); + + carDist.getUpper().setName("Caribbean High & Low"); + funcs.add(carDist.getUpper()); + chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 1f, carColor)); + + carDist.getLower().setName(null); + funcs.add(carDist.getLower()); + chars.add(chars.get(chars.size()-1)); + + PlotSpec plot = new PlotSpec(funcs, chars, " ", "Magnitude", "Incremental Rate (1/yr)"); + plot.setLegendInset(true); + + Range yRange = new Range(1e-6, 1e1); + + HeadlessGraphPanel gp = PlotUtils.initHeadless(); + + gp.drawGraphPanel(plot, false, true, xRange, yRange); + + PlotUtils.writePlots(outputDir, prefix, gp, 800, 750, true, true, false); + } private static double[] fractiles = {0d, 0.025, 0.16, 0.84, 0.975d, 1d}; private static String fractileLabel = "p[0,2.5,16,84,97.5,100]"; 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 223d1303..c6a22bda 100644 --- a/src/main/java/scratch/kevin/prvi25/figures/PRVI_SubductionSubSectPlots.java +++ b/src/main/java/scratch/kevin/prvi25/figures/PRVI_SubductionSubSectPlots.java @@ -34,6 +34,7 @@ import com.google.common.base.Preconditions; +import net.mahdilamb.colormap.Colors; import scratch.kevin.prvi25.FaultSystemLineIntegralCalculator; import static scratch.kevin.prvi25.figures.PRVI_Paths.*; @@ -46,7 +47,8 @@ public class PRVI_SubductionSubSectPlots { public static void main(String[] args) throws IOException { LogicTreeBranch branch = PRVI25_LogicTreeBranch.DEFAULT_SUBDUCTION_INTERFACE; - PRVI25_SubductionScalingRelationships scale = PRVI25_SubductionScalingRelationships.LOGA_C4p0; + PRVI25_SubductionScalingRelationships scale = PRVI25_SubductionScalingRelationships.AVERAGE; + String scaleLabel = "average scaling"; branch = branch.copy(); branch.setValue(scale); PRVI25_InvConfigFactory factory = new PRVI25_InvConfigFactory(); @@ -56,10 +58,10 @@ public static void main(String[] args) throws IOException { Preconditions.checkState(outputDir.exists() || outputDir.mkdir()); Font font = new Font(Font.SANS_SERIF, Font.BOLD, 24); - DecimalFormat magDF = new DecimalFormat("0.00"); + DecimalFormat magDF = new DecimalFormat("0.0"); CPT magCPT = GMT_CPT_Files.SEQUENTIAL_BATLOW_UNIFORM.instance().rescale(7.5d, 9d); - CPT slipCPT = GMT_CPT_Files.SEQUENTIAL_BATLOW_UNIFORM.instance().rescale(0d, 5d); + CPT slipCPT = GMT_CPT_Files.SEQUENTIAL_BATLOW_UNIFORM.instance().rescale(0d, 4d); CPT slipUncertCPT = GMT_CPT_Files.SEQUENTIAL_BATLOW_UNIFORM.instance().rescale(0d, 2d); CPT rakeCPT = GMT_CPT_Files.SEQUENTIAL_NAVIA_UNIFORM.instance().rescale(0d, 90d); @@ -101,14 +103,35 @@ public static void main(String[] args) throws IOException { double annX = xRange.getLowerBound() + 0.97*xRange.getLength(); double annY = yRange.getLowerBound() + 0.95*yRange.getLength(); - mapMaker.plotSectScalars(minMags, magCPT, "Minimum Magnitude ("+scale.getShortName()+")"); + mapMaker.plotSectScalars(minMags, magCPT, "Minimum Magnitude ("+scaleLabel+")"); XYTextAnnotation rangeAnn = new XYTextAnnotation("["+magDF.format(minMin)+", "+magDF.format(maxMin)+"]", annX, annY); rangeAnn.setTextAnchor(TextAnchor.TOP_RIGHT); rangeAnn.setFont(font); mapMaker.addAnnotation(rangeAnn); mapMaker.plot(outputDir, "subduction_min_mag_"+fm.getFilePrefix(), fmName); - mapMaker.plotSectScalars(maxMags, magCPT, "Maximum Magnitude ("+scale.getShortName()+")"); + Font subInterfaceFont = new Font(Font.SANS_SERIF, Font.BOLD, 18); + XYTextAnnotation hspAnn = new XYTextAnnotation("Northern Hispaniola", -69.8, 20.05); + hspAnn.setRotationAngle(mapMaker.getRotationAngleCorrectedForAspectRatio(Math.toRadians(18))); + hspAnn.setTextAnchor(TextAnchor.BASELINE_CENTER); + hspAnn.setFont(subInterfaceFont); + mapMaker.addAnnotation(hspAnn); + + XYTextAnnotation prviAnn = new XYTextAnnotation("Puerto Rico-Virgin Islands", -65.5, 20.05); +// prviAnn.setRotationAngle(mapMaker.getRotationAngleCorrectedForAspectRatio(Math.toRadians(20))); + prviAnn.setTextAnchor(TextAnchor.BASELINE_CENTER); + prviAnn.setFont(subInterfaceFont); + mapMaker.addAnnotation(prviAnn); + + XYTextAnnotation laAnn = new XYTextAnnotation("Lesser Antilles", -62, 19.55); + laAnn.setRotationAngle(mapMaker.getRotationAngleCorrectedForAspectRatio(Math.toRadians(23))); + laAnn.setTextAnchor(TextAnchor.BASELINE_CENTER); + laAnn.setFont(subInterfaceFont); + mapMaker.addAnnotation(laAnn); + + mapMaker.plot(outputDir, "subduction_min_mag_"+fm.getFilePrefix()+"_names", fmName); + + mapMaker.plotSectScalars(maxMags, magCPT, "Maximum Magnitude ("+scaleLabel+")"); mapMaker.clearAnnotations(); rangeAnn = new XYTextAnnotation("["+magDF.format(minMax)+", "+magDF.format(maxMax)+"]", annX, annY); rangeAnn.setTextAnchor(TextAnchor.TOP_RIGHT); @@ -116,6 +139,18 @@ public static void main(String[] args) throws IOException { mapMaker.addAnnotation(rangeAnn); mapMaker.plot(outputDir, "subduction_max_mag_"+fm.getFilePrefix(), fmName); + Font interfaceFont = new Font(Font.SANS_SERIF, Font.BOLD, 26); + XYTextAnnotation carAnn = new XYTextAnnotation("Caribbean Trench", -65.5, 20.05); + carAnn.setTextAnchor(TextAnchor.BASELINE_CENTER); + carAnn.setFont(interfaceFont); + mapMaker.addAnnotation(carAnn); + XYTextAnnotation mueAnn = new XYTextAnnotation("Muertos Trough", -68, 17.2); + mueAnn.setTextAnchor(TextAnchor.TOP_CENTER); + mueAnn.setFont(interfaceFont); + mapMaker.addAnnotation(mueAnn); + + mapMaker.plot(outputDir, "subduction_max_mag_"+fm.getFilePrefix()+"_names", fmName); + mapMaker.clearAnnotations(); for (PRVI25_SubductionDeformationModels dm : PRVI25_SubductionDeformationModels.values()) { @@ -142,10 +177,10 @@ public static void main(String[] args) throws IOException { 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); + rakeArrows.addAll(buildRakeArrows(sect, 5d)); + mapMaker.plotArrows(rakeArrows, 20d, Colors.tab_red.darker(), 2f); mapMaker.setFillArrowheads(true); - mapMaker.plotLines(rakeArrows, Color.BLACK, 2f); +// mapMaker.plotLines(rakeArrows, Color.BLACK, 2f); mapMaker.plotSectScalars(slips, slipCPT, dm.getShortName()+" Slip Rate (mm/yr)"); mapMaker.plot(outputDir, "subduction_slip_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), fmName+", "+dmName); diff --git a/src/main/java/scratch/kevin/prvi25/figures/ProxyFaultRepresentationFigures.java b/src/main/java/scratch/kevin/prvi25/figures/ProxyFaultRepresentationFigures.java index eb83c913..62b11d48 100644 --- a/src/main/java/scratch/kevin/prvi25/figures/ProxyFaultRepresentationFigures.java +++ b/src/main/java/scratch/kevin/prvi25/figures/ProxyFaultRepresentationFigures.java @@ -136,7 +136,7 @@ public static void main(String[] args) throws IOException { int rupID = -1; double maxRupMag = 0d; List allRups = new ArrayList<>(rupSet.getRupturesForParentSection(parentID)); -// Collections.reverse(allRups); // this makes it the left side of Anegada zones + Collections.reverse(allRups); // this makes it the left side of Anegada zones for (int rupIndex : allRups) { double mag = rupSet.getMagForRup(rupIndex); if ((float)mag < (float)maxRupMag) diff --git a/src/main/java/scratch/kevin/prvi25/figures/RupSetStatsTexWriter.java b/src/main/java/scratch/kevin/prvi25/figures/RupSetStatsTexWriter.java index f1a9f256..e58c9a0e 100644 --- a/src/main/java/scratch/kevin/prvi25/figures/RupSetStatsTexWriter.java +++ b/src/main/java/scratch/kevin/prvi25/figures/RupSetStatsTexWriter.java @@ -117,6 +117,8 @@ private static void write(FileWriter fw, String texPrefix, PRVI25_InvConfigFacto ClusterRuptures cRups = avgRupSet.requireModule(ClusterRuptures.class); BinaryRuptureProbabilityCalc exclusion = PRVI25_InvConfigFactory.getExclusionModel( avgRupSet, avgBranch, cRups); + double myMaxMin = 0d; + double myMinMin = Double.POSITIVE_INFINITY; for (int s=0; s linearAnns = circleFaultAnns(targetSubSects, "Main Ridge", solSlips, targetSlips, TextAnchor.TOP_LEFT); linearAnns.addAll(circleFaultAnns(targetSubSects, "Bouillante Montserrat", solSlips, targetSlips, TextAnchor.BASELINE_RIGHT)); + Range linearRange = new Range(0d, 13d); + Range logRange = new Range(1e-1, 2e1); + plotScatter(solOutputDir, prefix+"_scatter", targetSlips, solSlips, - targetAverageSlipsByParent, solAverageSlipsByParent, linearAnns, null); + targetAverageSlipsByParent, solAverageSlipsByParent, + linearRange, logRange, linearAnns, null); CSVFile slipsCSV = new CSVFile<>(true); @@ -301,6 +314,184 @@ private static void plotCrustal(File dmOutputDir, File solOutputDir, FaultSystem } } + private static void plotSubduction(File outputDir, FaultSystemSolution smallSol, FaultSystemSolution largeSol) throws IOException { + double weightLarge = PRVI25_SubductionFaultModels.PRVI_SUB_FM_LARGE.getNodeWeight(null); + double weightSmall = PRVI25_SubductionFaultModels.PRVI_SUB_FM_SMALL.getNodeWeight(null); + + Preconditions.checkState(smallSol.getRupSet().getNumSections() == largeSol.getRupSet().getNumSections()); + + int[] commonParents = { + FaultSectionUtils.findParentSectionID(smallSol.getRupSet().getFaultSectionDataList(), "Hispaniola"), + FaultSectionUtils.findParentSectionID(smallSol.getRupSet().getFaultSectionDataList(), "Muertos") + }; + + List combSects = new ArrayList<>(); + List combSectSolSlipRates = new ArrayList<>(); + List combSectTargetSlipRates = new ArrayList<>(); + Map combAverageSlipsByParent = new HashMap<>(); + Map combAverageTargetByParent = new HashMap<>(); + Map idToCombMap = new HashMap<>(); + + List absValues = new ArrayList<>(); + List scaledToRates = new ArrayList<>(); + + for (boolean small : new boolean[] {false,true}) { + List sects; + SolutionSlipRates solSlipsModule; + SectSlipRates targets; + double weight; + if (small) { + targets = smallSol.getRupSet().getSectSlipRates(); + sects = smallSol.getRupSet().getFaultSectionDataList(); + solSlipsModule = smallSol.requireModule(SolutionSlipRates.class); + weight = weightSmall/(weightSmall+weightLarge); + } else { + targets = largeSol.getRupSet().getSectSlipRates(); + sects = largeSol.getRupSet().getFaultSectionDataList(); + solSlipsModule = largeSol.requireModule(SolutionSlipRates.class); + weight = weightLarge/(weightSmall+weightLarge); + } + + List solSlips = new ArrayList<>(); + for (int i=0; i> solSlipsByParent = new HashMap<>(); + for (FaultSection sect : sects) { + int parentID = sect.getParentSectionId(); + if (!solSlipsByParent.containsKey(parentID)) + solSlipsByParent.put(parentID, new ArrayList<>()); + solSlipsByParent.get(parentID).add(solSlipsModule.get(sect.getSectionId())*1e3); + } + Map> targetSlipsByParent = new HashMap<>(); + for (int i=0; i()); + targetSlipsByParent.get(parentID).add(slip); + } + for (int s=0; s solAverageSlipsByParent = solSlipsByParent.entrySet().stream() + .collect(Collectors.toMap(Map.Entry::getKey, + entry -> entry.getValue().stream().mapToDouble(Double::doubleValue).average().orElse(0.0) + )); + Map targetAverageSlipsByParent = targetSlipsByParent.entrySet().stream() + .collect(Collectors.toMap(Map.Entry::getKey, + entry -> entry.getValue().stream().mapToDouble(Double::doubleValue).average().orElse(0.0) + )); + for (int parent : solAverageSlipsByParent.keySet()) { + double parentSolSlip = solAverageSlipsByParent.get(parent); + double parentTargetSlip = targetAverageSlipsByParent.get(parent); + if (Ints.contains(commonParents, parent)) { + double prev = combAverageSlipsByParent.containsKey(parent) ? combAverageSlipsByParent.get(parent) : 0d; + combAverageSlipsByParent.put(parent, prev + weight*parentSolSlip); + combAverageTargetByParent.put(parent, prev + weight*parentTargetSlip); + } else { + // unique fake parent ID + int fakeParent = -solAverageSlipsByParent.size(); + combAverageSlipsByParent.put(fakeParent, parentSolSlip); + combAverageTargetByParent.put(fakeParent, parentTargetSlip); + } + } + } + + GeographicMapMaker mapMaker = new GeographicMapMaker(PRVI_SubductionSubSectPlots.plotReg, combSects); + mapMaker.setWriteGeoJSON(false); + mapMaker.setFillSurfaces(true); + mapMaker.setSectTraceChar(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.DARK_GRAY)); + mapMaker.setSectOutlineChar(new PlotCurveCharacterstics(PlotLineType.SOLID, 1f, Color.GRAY)); + + 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); + + List sortables = new ArrayList<>(combSectSolSlipRates.size()); + for (int i=0; i slipsCSV = new CSVFile<>(true); + + slipsCSV.addLine("", "Absolute slip rate misfit (mm/yr)", "Misfit percentage of rate"); + DecimalFormat pDF = new DecimalFormat("0.#%"); + DecimalFormat twoDigits = new DecimalFormat("0.00"); + + FileWriter slipTEX = new FileWriter(new File(outputDir, "sub_slip_misfits.tex")); + + for (boolean max : new boolean[] {false,true}) { + String label = "Subsection"; + if (max) + label += ", maximum"; + else + label += ", average"; + List line = new ArrayList<>(); + line.add(label); + double abs, scaledToRate; + String texPrefix = "SubductionSubsect"; + if (max) { + texPrefix += "Max"; + abs = absValues.stream().mapToDouble(D->D).max().getAsDouble(); + scaledToRate = scaledToRates.stream().mapToDouble(D->D).max().getAsDouble(); + } else { + texPrefix += "Avg"; + abs = absValues.stream().mapToDouble(D->D).average().getAsDouble(); + scaledToRate = scaledToRates.stream().mapToDouble(D->D).average().getAsDouble(); + } + texPrefix += "SlipMsft"; + slipTEX.write(LaTeXUtils.defineValueCommand(texPrefix, twoDigits.format(abs))+"\n"); + slipTEX.write(LaTeXUtils.defineValueCommand(texPrefix+"Pct", pDF.format(scaledToRate))+"\n"); + line.add(twoDigits.format(abs)); + line.add(pDF.format(scaledToRate)); + slipsCSV.addLine(line); + } + + slipTEX.close(); + slipsCSV.writeToFile(new File(outputDir, "sub_slip_misfits.csv")); + } + private static List circleFaultAnns(List sects, String name, List solSlips, List targetSlips, TextAnchor anchor) { List anns = new ArrayList<>(); @@ -392,19 +583,18 @@ private static List percentDifference(List target, List return ret; } private static void plotScatter(File outputDir, String prefix, List targetSlips, List solSlips, - Map targetParentSlips, Map solParentSlips) throws IOException { - plotScatter(outputDir, prefix, targetSlips, solSlips, targetParentSlips, solParentSlips, null, null); + Map targetParentSlips, Map solParentSlips, + Range linearRange, Range logRange) throws IOException { + plotScatter(outputDir, prefix, targetSlips, solSlips, targetParentSlips, solParentSlips, linearRange, logRange, null, null); } private static void plotScatter(File outputDir, String prefix, List targetSlips, List solSlips, Map targetParentSlips, Map solParentSlips, + Range linearRange, Range logRange, List linearAnns, List logAnns) 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()); diff --git a/src/main/java/scratch/kevin/ucerf3/PureScratch.java b/src/main/java/scratch/kevin/ucerf3/PureScratch.java index d1ec73cb..5ec0c7ff 100644 --- a/src/main/java/scratch/kevin/ucerf3/PureScratch.java +++ b/src/main/java/scratch/kevin/ucerf3/PureScratch.java @@ -3073,13 +3073,48 @@ private static void test336() throws IOException { System.out.println("Slip: "+(float)slip+" -> "+(float)projected); } + private static void test337() throws IOException { + int num = 100000000; + int numVals = 1000; + double[] randBases = new double[numVals]; + double[] randExponents = new double[numVals]; + for (int i=0; i