Skip to content

Commit

Permalink
initial interpolation scheme
Browse files Browse the repository at this point in the history
  • Loading branch information
Kevin Milner committed Jun 3, 2024
1 parent d728ca4 commit 08aedbb
Show file tree
Hide file tree
Showing 2 changed files with 86 additions and 3 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@

import org.jfree.chart.annotations.XYTextAnnotation;
import org.jfree.chart.ui.TextAnchor;
import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc;
import org.opensha.commons.geo.Location;
import org.opensha.commons.geo.LocationList;
import org.opensha.commons.geo.LocationUtils;
Expand Down Expand Up @@ -62,6 +63,8 @@ public static void main(String[] args) throws IOException {
boolean[] fullRateBools = { true, false };
boolean[] largePolyBools = { true, false };

boolean interpolate = true;

File debugDir = new File("/tmp/sub_fm_dm_debug");
Preconditions.checkState(debugDir.exists() || debugDir.mkdir());

Expand Down Expand Up @@ -158,6 +161,12 @@ public static void main(String[] args) throws IOException {

int prevID = -1;

List<Double> minisectionDASs = new ArrayList<>();
ArbitrarilyDiscretizedFunc dasSlipFunc = new ArbitrarilyDiscretizedFunc();
ArbitrarilyDiscretizedFunc dasSlipUncertFunc = new ArbitrarilyDiscretizedFunc();
ArbitrarilyDiscretizedFunc dasRakeFunc = new ArbitrarilyDiscretizedFunc();
double runningDAS = 0d;

for (Feature feature : features) {
Preconditions.checkState(feature.geometry.type == GeoJSON_Type.MultiLineString);
MultiLineString geometry = (MultiLineString)feature.geometry;
Expand Down Expand Up @@ -226,16 +235,43 @@ public static void main(String[] args) throws IOException {
double slipUncert = props.getDouble(rateUncertPropName, Double.NaN);
double rake = props.getDouble(rakePropName, Double.NaN);

// TODO convert to plane rates
if (runningDAS == 0d) {
dasSlipFunc.set(0d, slip);
dasSlipUncertFunc.set(0d, slipUncert);
dasRakeFunc.set(0d, rake);
}

double myStartDAS = runningDAS;

for (int i=1; i<upper.size(); i++) {
Location startLoc = upper.get(i-1);
Location endLoc = upper.get(i);
int minisectionID = minisectionRecs.size();
minisectionRecs.add(new MinisectionSlipRecord(outID, minisectionID, startLoc, endLoc, rake, slip, slipUncert));
double len = LocationUtils.horzDistanceFast(startLoc, endLoc);
minisectionDASs.add(runningDAS + 0.5*len); // middle
runningDAS += len;
}

double myEndDAS = runningDAS;
double middleDAS = myStartDAS + 0.5*(myEndDAS - myStartDAS);

System.out.println("\tDAS for "+sectName+": "+middleDAS);

dasSlipFunc.set(middleDAS, slip);
dasSlipUncertFunc.set(middleDAS, slipUncert);
dasRakeFunc.set(middleDAS, rake);

prevID = sectID;
}

System.out.println("Total DAS for "+name+": "+runningDAS);

// add end
dasSlipFunc.set(runningDAS, dasSlipFunc.getY(dasSlipFunc.size()-1));
dasSlipUncertFunc.set(runningDAS, dasSlipUncertFunc.getY(dasSlipUncertFunc.size()-1));
dasRakeFunc.set(runningDAS, dasRakeFunc.getY(dasRakeFunc.size()-1));

// build stitched feature
MultiLineString geometry = new MultiLineString(List.of(stitchedUpperTrace, stitchedLowerTrace));
FeatureProperties props = new FeatureProperties();
Expand All @@ -244,6 +280,21 @@ public static void main(String[] args) throws IOException {
props.set("PrimState", "PR");
outFeatures.add(new Feature(outID, geometry, props));

if (interpolate) {
// interpolate slips and rakes
Preconditions.checkState(minisectionDASs.size() == minisectionRecs.size());
for (int i=0; i<minisectionRecs.size(); i++) {
MinisectionSlipRecord rec = minisectionRecs.get(i);
double das = minisectionDASs.get(i);
Preconditions.checkState(das >= 0 && das <= runningDAS);
minisectionRecs.set(i, new MinisectionSlipRecord(rec.parentID, rec.minisectionID,
rec.startLoc, rec.endLoc,
dasRakeFunc.getInterpolatedY(das),
dasSlipFunc.getInterpolatedY(das),
dasSlipUncertFunc.getInterpolatedY(das)));
}
}

minisectionRecsMap.put(outID, minisectionRecs);
}

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
import org.opensha.sha.earthquake.rupForecastImpl.prvi25.logicTree.PRVI25_SubductionScalingRelationships;
import org.opensha.sha.faultSurface.FaultSection;

public class SubductionMinMaxMagPlots {
public class PRVI_SubductionSubSectPlots {

public static void main(String[] args) throws IOException {
LogicTreeBranch<LogicTreeNode> branch = PRVI25_LogicTreeBranch.DEFAULT_SUBDUCTION_INTERFACE;
Expand All @@ -39,7 +39,12 @@ public static void main(String[] args) throws IOException {
DecimalFormat magDF = new DecimalFormat("0.00");

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 slipUncertCPT = GMT_CPT_Files.SEQUENTIAL_BATLOW_UNIFORM.instance().rescale(0d, 2d);
CPT rakeCPT = GMT_CPT_Files.SEQUENTIAL_NAVIA_UNIFORM.instance().rescale(0d, 90d);
for (PRVI25_SubductionFaultModels fm : PRVI25_SubductionFaultModels.values()) {
if (fm.getNodeWeight(branch) == 0d)
continue;
branch = branch.copy();
branch.setValue(fm);
branch.requireValue(PRVI25_SubductionDeformationModels.class).build(fm);
Expand Down Expand Up @@ -74,14 +79,41 @@ public static void main(String[] args) throws IOException {
mapMaker.addAnnotation(rangeAnn);
mapMaker.plot(outputDir, "subduction_min_mag_"+fm.getFilePrefix(), " ");


mapMaker.plotSectScalars(maxMags, magCPT, "Maximum Magnitude ("+scale.getShortName()+")");
mapMaker.clearAnnotations();
rangeAnn = new XYTextAnnotation("["+magDF.format(minMax)+", "+magDF.format(maxMax)+"]", annX, annY);
rangeAnn.setTextAnchor(TextAnchor.TOP_RIGHT);
rangeAnn.setFont(font);
mapMaker.addAnnotation(rangeAnn);
mapMaker.plot(outputDir, "subduction_max_mag_"+fm.getFilePrefix(), " ");

mapMaker.clearAnnotations();

for (PRVI25_SubductionDeformationModels dm : PRVI25_SubductionDeformationModels.values()) {
if (dm.getNodeWeight(branch) == 0d)
continue;
sects = dm.build(fm);

mapMaker.setFaultSections(sects);
List<Double> slips = new ArrayList<>(sects.size());
List<Double> slipUncerts = new ArrayList<>(sects.size());
List<Double> rakes = new ArrayList<>(sects.size());
for (FaultSection sect : sects) {
slips.add(sect.getOrigAveSlipRate());
slipUncerts.add(sect.getOrigSlipRateStdDev());
rakes.add(sect.getAveRake());
}

mapMaker.plotSectScalars(slips, slipCPT, dm.getShortName()+" Slip Rate (mm/yr)");
mapMaker.plot(outputDir, "subduction_slip_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), " ");

mapMaker.plotSectScalars(slipUncerts, slipUncertCPT, dm.getShortName()+" Slip Rate Uncertainty (mm/yr)");
mapMaker.plot(outputDir, "subduction_slip_uncert_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), " ");

mapMaker.plotSectScalars(rakes, rakeCPT, dm.getShortName()+" Rake");
mapMaker.plot(outputDir, "subduction_rake_"+fm.getFilePrefix()+"_"+dm.getFilePrefix(), " ");

}
}
}

Expand Down

0 comments on commit 08aedbb

Please sign in to comment.