add isochrone script #896
add isochrone script #896
Conversation
This script computes the isochrones (on-fault contours that are associated with acceleration waveforms) of a certain surface receiver.
|
I would recommend that we follow the standard Python naming/style conventions: https://peps.python.org/pep-0008/? |
|
I suppose that means applying a linter (e.g. black but you probably already did that) and checking style errors detected by flake8? |
|
The most obvious issue is the naming scheme: In Python, functions always use snake_case instead of CamelCase. |
| help='provide path+filename-fault.xdmf; only needed when the path differs from the -surface.xdmf file') | ||
|
|
||
| parser.add_argument('--output', choices=["numpy","xdmf","both"], default="xdmf", | ||
| help='choose the output format') |
There was a problem hiding this comment.
For consistency with the rest: help='output format'
| parser.add_argument('--output', choices=["numpy","xdmf","both"], default="xdmf", | ||
| help='choose the output format') | ||
|
|
||
| parser.add_argument('--sliprateThreshold', nargs=1, default=([0.05]), metavar=('peak sliprate threshold in m/s'), |
There was a problem hiding this comment.
--slipRateThreshold ?
| help='choose the output format') | ||
|
|
||
| parser.add_argument('--sliprateThreshold', nargs=1, default=([0.05]), metavar=('peak sliprate threshold in m/s'), | ||
| help='fault cells below the treshold are assumed to not radiate waves' , type=float) |
| parser.add_argument('--sliprateThreshold', nargs=1, default=([0.05]), metavar=('peak sliprate threshold in m/s'), | ||
| help='fault cells below the treshold are assumed to not radiate waves' , type=float) | ||
|
|
||
| parser.add_argument('--sliprateParameters', choices=["absolute","components","both"], default="both", |
There was a problem hiding this comment.
--slipRateParameters?
| help="absolute: compute isochrones from peak times of absolute slip rates;"+ | ||
| " components: peak slip rate times are computed for SRs and SRd separately") | ||
|
|
||
| parser.add_argument('--medianFilterRadius', nargs=1, default=([0.]), metavar=('radius of the median filter in m'), |
There was a problem hiding this comment.
maybe remove the default, and add required=True, forcing the user to have this (maybe) important decision?
|
|
||
| def ApproximateHypocenter(absoluteSliprates, faultCells, slipRateThreshold=args.sliprateThreshold[0]): | ||
|
|
||
| slippingFault = np.where(absoluteSliprates > slipRateThreshold, 1, 0) |
There was a problem hiding this comment.
wouldn't it be more simple to use the min of RT variable?
There was a problem hiding this comment.
yes, I rewrote the function
| faultXdmf = seissolxdmf.seissolxdmf(args.faultXdmf) | ||
| faultGeom = faultXdmf.ReadGeometry() | ||
| faultConnect = faultXdmf.ReadConnect() | ||
| faultCells = ComputeTriangleMidpoints(faultGeom, faultConnect) |
There was a problem hiding this comment.
faultMidPoints = ...
| print(f"sliprates.shape: {sliprates.shape}") | ||
|
|
||
| stop1 = timeit.default_timer() | ||
| print(f"Time to load data: {np.round(stop1 - start, 2)}") |
There was a problem hiding this comment.
print(f"Time to load data: {stop1 - start:.2f}"
| print(f"Hypocenter approximated at: {np.round(hypocenter)}") | ||
|
|
||
| avgSwaveVelocity = np.linalg.norm(receiverCoords - hypocenter) / receiverPwaveTraveltime / np.sqrt(3) # Assumes a Poisson solid | ||
| print(f"Average S-wave velocity between fault and receiver: {np.round(avgSwaveVelocity,2)}") |
There was a problem hiding this comment.
between hypocenter and receiver
| stop3 = timeit.default_timer() | ||
| print("Saving output...") | ||
|
|
||
| prefix = f"isochrones_{np.int(receiverCoords[0])}_{np.int(receiverCoords[1])}_{np.int(receiverCoords[2])}" |
There was a problem hiding this comment.
f"isochrones_{receiverCoords[0]:.0f}_(...)"
| help='provide path+filename-fault.xdmf; only needed when the path differs from the -surface.xdmf file') | ||
|
|
||
| parser.add_argument('--output', choices=["numpy","xdmf","both"], default="xdmf", | ||
| help='choose the output format') |
| parser.add_argument('--output', choices=["numpy","xdmf","both"], default="xdmf", | ||
| help='choose the output format') | ||
|
|
||
| parser.add_argument('--sliprateThreshold', nargs=1, default=([0.05]), metavar=('peak sliprate threshold in m/s'), |
| help='choose the output format') | ||
|
|
||
| parser.add_argument('--sliprateThreshold', nargs=1, default=([0.05]), metavar=('peak sliprate threshold in m/s'), | ||
| help='fault cells below the treshold are assumed to not radiate waves' , type=float) |
| parser.add_argument('--sliprateThreshold', nargs=1, default=([0.05]), metavar=('peak sliprate threshold in m/s'), | ||
| help='fault cells below the treshold are assumed to not radiate waves' , type=float) | ||
|
|
||
| parser.add_argument('--sliprateParameters', choices=["absolute","components","both"], default="both", |
| help="absolute: compute isochrones from peak times of absolute slip rates;"+ | ||
| " components: peak slip rate times are computed for SRs and SRd separately") | ||
|
|
||
| parser.add_argument('--medianFilterRadius', nargs=1, default=([0.]), metavar=('radius of the median filter in m'), |
| faultXdmf = seissolxdmf.seissolxdmf(args.faultXdmf) | ||
| faultGeom = faultXdmf.ReadGeometry() | ||
| faultConnect = faultXdmf.ReadConnect() | ||
| faultCells = ComputeTriangleMidpoints(faultGeom, faultConnect) |
| ruptureTimes = faultXdmf.ReadData("RT", timeIndicesFault[1]-1) | ||
|
|
||
| if args.events[0] == 2: | ||
| ruptureTimes = ruptureTimes - faultXdmf.ReadData("RT", timeIndicesFault[0]-1) |
There was a problem hiding this comment.
this sets all RTs to zero that slipped before timeIndicesFault[0]
| print(f"sliprates.shape: {sliprates.shape}") | ||
|
|
||
| stop1 = timeit.default_timer() | ||
| print(f"Time to load data: {np.round(stop1 - start, 2)}") |
| print(f"Hypocenter approximated at: {np.round(hypocenter)}") | ||
|
|
||
| avgSwaveVelocity = np.linalg.norm(receiverCoords - hypocenter) / receiverPwaveTraveltime / np.sqrt(3) # Assumes a Poisson solid | ||
| print(f"Average S-wave velocity between fault and receiver: {np.round(avgSwaveVelocity,2)}") |
| stop3 = timeit.default_timer() | ||
| print("Saving output...") | ||
|
|
||
| prefix = f"isochrones_{np.int(receiverCoords[0])}_{np.int(receiverCoords[1])}_{np.int(receiverCoords[2])}" |
|
@krenzland I will follow the style guide in the future. Do you think it's necessary to update the function/variable names of this script? |
|
|
||
| def ComputeTriangleMidpoints(geom, connect): | ||
| """Generates an array with coordinates of triangle midpoints (same dimension as connect)""" | ||
| xyz = np.zeros_like(connect, dtype=float) |
There was a problem hiding this comment.
this line is not required
|
|
||
|
|
||
| def LoadSliprates(sliprateParameters): | ||
| sliprates = [] |
There was a problem hiding this comment.
sliprates = {}
sliprates["SRs"] = tXdmf.ReadData("SRs")[timeIndicesFault[0]:timeIndicesFault[1]].T
sliprates["SRd"] = tXdmf.ReadData("SRd")[timeIndicesFault[0]:timeIndicesFault[1]].T
etc.
I mean it would be a dictionary, with keys SRs, SRd, and ASR pointing to numpy array.
| faultConnect = faultXdmf.ReadConnect() | ||
| faultMidPoints = ComputeTriangleMidpoints(faultGeom, faultConnect) | ||
| dtFault = faultXdmf.ReadTimeStep() | ||
| timeIndicesFault = np.array(args.timeStepRange) |
There was a problem hiding this comment.
What about naming this array faultSlicingRange? indexBoundsFault?
|
|
||
| print("Loading data...") | ||
| surfaceXdmf = seissolxdmf.seissolxdmf(args.filename) | ||
| surfaceCells = ComputeTriangleMidpoints(surfaceXdmf.ReadGeometry(), surfaceXdmf.ReadConnect()) |
There was a problem hiding this comment.
rename to surfaceCellBarycenters? or freeSurfaceBarycenters?
| def PTraveltime(receiverIndex, trigger=0.01): | ||
| """Picks the P-wave arrival. This function is only suited for synthetic data with negligible noise""" | ||
|
|
||
| timeIndices = np.int_(timeIndicesFault * (dtFault / dtSurface)) |
There was a problem hiding this comment.
freeSurfaceSlicingRange?
| surfaceXdmf.ReadDataChunk(surfaceVariables[2], receiverIndex, 1)**2)[timeIndices[0]:timeIndices[1]] | ||
|
|
||
| threshold = np.amax(absWaveform)*trigger # trigger dictates the portion of the maximum ground velocity, | ||
| t = np.argmax(absWaveform>threshold)*dtSurface # which needs to be reached to pick the p-phase |
There was a problem hiding this comment.
return np.argmax(absWaveform>threshold)*dtSurface
| stop2 = timeit.default_timer() | ||
| print(f"Time to compute values: {stop2 - stop1:.2f}") | ||
|
|
||
| if args.medianFilterRadius[0] > 0.: |
There was a problem hiding this comment.
or
if args.medianFilterRadius:
|
|
||
| def MedianSmoothingParallel(parameter, xyz, radius=250., nprocs=4): | ||
|
|
||
| chunks = [tuple((parameter, xyz, i, radius)) for i in np.array_split(np.arange(parameter.shape[0]), nprocs)] |
There was a problem hiding this comment.
chunks = [(parameter, xyz, i, radius) for i in np.array_split(np.arange(parameter.shape[0]), nprocs)]
should be sufficient.
| return parameter_new | ||
|
|
||
|
|
||
| def MedianSmoothingParallel(parameter, xyz, radius=250., nprocs=4): |
There was a problem hiding this comment.
would be clear if parameter was called
peakSliprateTimes (or a slightly different name if you want to avoid the name of the global value)
| parameter = input_tuple[0] | ||
| xyz = input_tuple[1] | ||
| indices = input_tuple[2] | ||
| radius = input_tuple[3] |
There was a problem hiding this comment.
parameter, xyz, indices, radius = input_tuple
| radius = input_tuple[3] | ||
| parameter_new = np.zeros(indices.shape[0]) | ||
| for i, j in enumerate(indices): | ||
| parameter_new[i] = np.median(parameter[(xyz[:,0] < xyz[j,0]+radius) & (xyz[:,0] > xyz[j,0]-radius) & |
There was a problem hiding this comment.
for i, j in enumerate(indices):
xyz0= xyz[j,:]
distances = np.linalg.norm(xyz - xyz0, axis=1)
indices_within_sphere = np.where(distances <= radius)[0]
parameter_new[i] =np.median(parameters[indices_within_sphere])
There was a problem hiding this comment.
Yes, much cleaner
|
|
||
|
|
||
| def LoadSliprates(sliprateParameters): | ||
| sliprates = [] |
There was a problem hiding this comment.
But later array operations are performed on sliprates. If it was a dictionary, there would be a loop necessary each time.
E.g., here:
sliprates = np.where(sliprates >= args.slipRateThreshold[0], sliprates, 0.) # Set slip rates below the threshold to zero
peakSliprateTimes = np.argmax(sliprates, axis=2) * dtFault
| def PTraveltime(receiverIndex, trigger=0.01): | ||
| """Picks the P-wave arrival. This function is only suited for synthetic data with negligible noise""" | ||
|
|
||
| timeIndices = np.int_(timeIndicesFault * (dtFault / dtSurface)) |
| faultConnect = faultXdmf.ReadConnect() | ||
| faultMidPoints = ComputeTriangleMidpoints(faultGeom, faultConnect) | ||
| dtFault = faultXdmf.ReadTimeStep() | ||
| timeIndicesFault = np.array(args.timeStepRange) |
| stop2 = timeit.default_timer() | ||
| print(f"Time to compute values: {stop2 - stop1:.2f}") | ||
|
|
||
| if args.medianFilterRadius[0] > 0.: |
|
|
||
| def ComputeTriangleMidpoints(geom, connect): | ||
| """Generates an array with coordinates of triangle midpoints (same dimension as connect)""" | ||
| xyz = np.zeros_like(connect, dtype=float) |
| parameter = input_tuple[0] | ||
| xyz = input_tuple[1] | ||
| indices = input_tuple[2] | ||
| radius = input_tuple[3] |
| radius = input_tuple[3] | ||
| parameter_new = np.zeros(indices.shape[0]) | ||
| for i, j in enumerate(indices): | ||
| parameter_new[i] = np.median(parameter[(xyz[:,0] < xyz[j,0]+radius) & (xyz[:,0] > xyz[j,0]-radius) & |
There was a problem hiding this comment.
Yes, much cleaner
| return parameter_new | ||
|
|
||
|
|
||
| def MedianSmoothingParallel(parameter, xyz, radius=250., nprocs=4): |
|
|
||
| def MedianSmoothingParallel(parameter, xyz, radius=250., nprocs=4): | ||
|
|
||
| chunks = [tuple((parameter, xyz, i, radius)) for i in np.array_split(np.arange(parameter.shape[0]), nprocs)] |
|
|
||
| print("Loading data...") | ||
| surfaceXdmf = seissolxdmf.seissolxdmf(args.filename) | ||
| surfaceCells = ComputeTriangleMidpoints(surfaceXdmf.ReadGeometry(), surfaceXdmf.ReadConnect()) |
Codecov Report
❗ Your organization needs to install the Codecov GitHub app to enable full functionality. @@ Coverage Diff @@
## master #896 +/- ##
=======================================
Coverage 14.35% 14.35%
=======================================
Files 253 253
Lines 14253 14253
=======================================
Hits 2046 2046
Misses 12207 12207 📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
| stop2 = timeit.default_timer() | ||
| print(f"Time to compute values: {stop2 - stop1:.2f}") | ||
|
|
||
| if args.medianFilterRadius[0]: |
There was a problem hiding this comment.
does this really works?
I would have think only if args.medianFilterRadius: works.
There was a problem hiding this comment.
Yes, it works because it's a required parameter
There was a problem hiding this comment.
then it is always true.
There was a problem hiding this comment.
No, you deactivate it by setting it to zero
This script computes the isochrones (on-fault contours that are associated with acceleration waveforms) of a certain surface receiver.