BEAT: Added option to build GFs with EPcrust model #87
Conversation
mgovorcin
changed the title
There was a problem hiding this comment.
Thanks a lot for your very valuable contribution! This is greatly appreciated!
We still need a test and documentation for the shapefile parser, e.g. through an example shapefile and a short explanation how this could be easily created. This could be very effectively covered using the same example file. Not sure how messy it would become if we would include the shape-file parser function to beat import? Maybe its better to keep as stand alone? Whats your opinion?
These notes are also for my own reference and please do not feel pressured to fulfill my suggestions and comments within any time-frame. Again everything is greatly appreciated!
| import os, sys | ||
| import argparse | ||
|
|
||
| class faultSegments: |
There was a problem hiding this comment.
Awesome script! Thanks a lot!
Please change name formatting to pyrocko/BEAT conventions.
class FaultSegments:
...
def add_points_lonlat(self):
...
def import_fault_shape():
| def mfaultToBeat(faultSegments): | ||
| # default lower upper testvalue | ||
| priors = ['depth', 'dip','duration', 'east_shift', 'length', 'north_shift', 'nucleation_x', 'nucleation_y', 'rake', 'slip', 'strike', 'time', 'width'] | ||
| default_bounds = { |
There was a problem hiding this comment.
could you use the defaults in beat.config for consistency? or is there anything preventing that?
Ah I see you need to sort between extracted from shape and these bounds here ...
| for i in range(len(fault)): | ||
| faultXY = orthodrome.latlon_to_ne_numpy(fault[i][:,1],fault[i][:,0],event[1],event[0]) | ||
|
|
||
| fSegmentPointsLonLat = np.zeros((len(faultXY[0])-1,4)) |
There was a problem hiding this comment.
faultXY is a numpy array, so instead of len(faultXY[0]), its better to use:
faultXY.shape[0]
why the -1? Because of the closed polygons thus need to delete last line?
define nsegment_points = faultXY.shape[0] and then reuse instead of reevaluating many times.
| plt.title(" Fault segements to model in BEAT") | ||
|
|
||
|
|
||
| def import_faultShp(daShapefile): |
There was a problem hiding this comment.
need documentation somewhere how the shapefile needs to be structured to be readable with this script.
| else: | ||
| profile = ep_crust.get_profile(location.lat, location.lon) | ||
| logger.info('EPcrust profile on the coordinates: %.2f, %.2f' % (profile._lat, profile._lon)) | ||
| if gfc.replace_water: |
There was a problem hiding this comment.
can we reuse that functionality from the crust2 stuff? Thus could we create a function:
def remove_water from profile():
and use that in the crust2 and the ep_crust?
| tabulated_phases.append(gf.TPDef( | ||
| id='any_S', | ||
| definition='s,S,s\\,S\\')) | ||
|
|
||
| # surface waves | ||
| if 'slowest' in waveforms: |
There was a problem hiding this comment.
always adding slowest may become very expensive for teleseismic setups, wouldnt want to hardcode add it.
Why do you need it here?
| @@ -1840,20 +1885,38 @@ def choose_backend( | |||
|
|
|||
| pids = ['stored:' + wave for wave in waveforms] | |||
|
|
|||
| ################################################################################# | |||
There was a problem hiding this comment.
I remember you told me sth about that in our chats. What was the issue? Maybe we should expose these offset-factors to the gf_config somehow? If we cant figure it out we can also split the PR into that of your shape-file parser and that to the gf_config rework.
| @@ -22,6 +22,7 @@ | |||
| from pyrocko.guts_array import Array | |||
|
|
|||
| from pyrocko import crust2x2, gf, cake, orthodrome, trace, util | |||
| from pyrocko.dataset import ep_crust | |||
There was a problem hiding this comment.
need a try-except surrounding that in case the pyrocko version does not include it, and then set the ep_crust to False or better raise an error with a message that pyrocko needs updating ... Its still not merged there right? Will also look through your PR there, .... soonish.
add the following lines below the "gf_config: !beat.SeismicGFConfig" in your config.yaml file
use_crust2: false #need to deselect this option as you want to use epcrust
use_epcrust: true
epcrust_radius: 20 # added this option to get the average EPcrust profile in the defined radius due to the high resolution of the model (0.5°x0.5°), use 0 if you want to get the EPcrust profile close to your reference location
In order to work pyrocko version needs to have the following:
modified cake.py to accept ep_crust and
ep_crust.py, EPcrust_0_5.txt, Ice_thickness_0_5.txt in the folder pyrocko/dataset
pyrocko.zip
######################################################
Added rough script to import multisegment fault from shapefile and export configuration file for "priors" - Maybe you can see how to get pieces of this code and include it somewhere else in beat, such as heat