diff --git a/beat/apps/importFaultShp.py b/beat/apps/importFaultShp.py new file mode 100644 index 00000000..2b7118a0 --- /dev/null +++ b/beat/apps/importFaultShp.py @@ -0,0 +1,239 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Jun 1 18:22:36 2021 + +@author: Marin Govorcin +""" + +import numpy as np +import os, sys +import argparse + +class faultSegments: + def __init__(self, segmentPointsLonLat= None, segmentPointsXY= None, segmentMidPoints= None, segmentStrike= None, segmentLength= None): + if segmentPointsLonLat is None: + self.segmentPointsLonLat = [] + else: + self.segmentPointsLonLat = segmentPointsLonLat + if segmentPointsXY is None: + self.segmentPointsXY = [] + else: + self.segmentPointsXY = segmentPointsXY + if segmentMidPoints is None: + self.segmentMidPoints = [] + else: + self.segmentMidPoints = segmentMidPoints + if segmentStrike is None: + self.segmentStrike = [] + else: + self.segmentStrike = segmentStrike + if segmentLength is None: + self.segmentLength = [] + else: + self.segmentLength = segmentLength + + def addPointsLonLat(self,data): + self.segmentPointsLonLat.append(data) + + def addPointsXY(self,data): + self.segmentPointsXY.append(data) + + def addMidPoints(self,data): + self.segmentMidPoints.append(data) + + def addStrikes(self,data): + self.segmentStrike.append(data) + + def addLengths(self,data): + self.segmentLength.append(data) + + def plotFault(self): + from matplotlib import pyplot as plt + + #plot fault segments + z=1 + for i in range(len(self.segmentPointsXY)): + for j in range(len(self.segmentPointsXY[i])): + plt.plot(self.segmentPointsXY[i][j,:].reshape(2,2)[0],self.segmentPointsXY[i][j,:].reshape(2,2)[1],'-k') + plt.plot(self.segmentMidPoints[i][j,0],self.segmentMidPoints[i][j,1],'bo') + plt.text(self.segmentMidPoints[i][j,0]+100,self.segmentMidPoints[i][j,1]+100, 'F{} , Strike: '.format(z) + str(int(self.segmentStrike[i][j]))+'°') + z = z +1 + #plot event location + plt.plot(0,0,'r+', label = 'Event location',ms=14) + plt.xlabel("East shift [m]") + plt.ylabel("North shift [m]") + plt.title(" Fault segements to model in BEAT") + + +def import_faultShp(daShapefile): + try: + from osgeo import ogr + print('Import of ogr from osgeo worked!\n') + except: + print('Import of ogr from osgeo failed\n\n') + + # import data from shapefile + driver = ogr.GetDriverByName('ESRI Shapefile') + data = driver.Open(daShapefile,0) + + # Check to see if shapefile is found. + if data is None: + print('Could not open %s' % daShapefile) + else: + print('Opened %s' % daShapefile) + layer = data.GetLayer() + featureCount = layer.GetFeatureCount() + print("Number of features in %s: %d" % (os.path.basename(daShapefile),featureCount)) + + #extract Fault Point Coordinates from shapefile + faultPoints = {} + + for i in range(featureCount): + feature = layer.GetFeature(i) + featureGeometry = feature.GetGeometryRef() + fnPoints = featureGeometry.GetPointCount() + featurePoints = np.zeros((fnPoints,2)) + + #extract points from shapefile features + for j in range(fnPoints): + featurePoints[j,:] = np.asarray(featureGeometry.GetPoint(j)[0:2]) + + faultPoints[i] = featurePoints + + # Output array of fault points + return faultPoints + + +def getFaultSegments(event, faultSource): + from pyrocko import orthodrome + #Get fault points coordinates in local coordinate system with respect to the event location + + # Get point coordinates from shapefile + # event :: list [lon, lat] + # faultSource :: str "filename.shp" + fault = import_faultShp(faultSource) + + fSegments = faultSegments() + + #get Points from fault features - lines + 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)) + fSegmentPointsXY = np.zeros((len(faultXY[0])-1,4)) + fSegmentMidPoint = np.zeros((len(faultXY[0])-1,2)) + fSegmentStrike = np.zeros((len(faultXY[0])-1,1)) + fSegmentLength = np.zeros((len(faultXY[0])-1,1)) + + + for j in range(len(faultXY[0])-1): + fSegmentPointsLonLat[j,:] = [fault[i][j][0], fault[i][j+1][0], fault[i][j][1], fault[i][j+1][1]] + fSegmentPointsXY[j,:] = [faultXY[1][j]*-1, faultXY[1][j+1]*-1, faultXY[0][j]*-1, faultXY[0][j+1]*-1] # x1 x2 y1 y2 + fSegmentMidPoint[j,:] = [(faultXY[1][j] + faultXY[1][j+1])/-2, (faultXY[0][j] + faultXY[0][j+1])/-2] + fSegmentStrike[j,:] = orthodrome.azimuth(fault[i][j,1],fault[i][j,0],fault[i][j+1,1],fault[i][j+1,0]) + fSegmentLength[j,:] = orthodrome.distance_accurate15nm(fault[i][j,1],fault[i][j,0],fault[i][j+1,1],fault[i][j+1,0]) + + fSegments.addPointsLonLat(fSegmentPointsLonLat) + fSegments.addPointsXY(fSegmentPointsXY) + fSegments.addMidPoints(fSegmentMidPoint) #[East_shift, North_shift] + fSegments.addStrikes(fSegmentStrike) + fSegments.addLengths(fSegmentLength) + + return fSegments + +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 = { + 'rake' : [0, 180, 90], + 'dip' : [0, 90, 45], + 'depth' : [0, 30, 15], + 'slip' : [0, 5, 2.5], + 'width' : [0, 10, 5], + 'duration' : [0, 30, 15], + 'nucleation_x' : [-1, 1, 0], + 'nucleation_y' : [-1, 1, 0], + 'time' : [-5, 5, 0]} + + file=open('beatPriors_multiSegment.txt','w') + + midPoints = np.vstack(faultSegments.segmentMidPoints) + east_shift = midPoints[:,0] / 1000 # in km + north_shift = midPoints[:,1] / 1000 # in km + strike = np.vstack(faultSegments.segmentStrike) + length = np.vstack(faultSegments.segmentLength) / 1000 # in km + + bounds = {'east_shift':east_shift, 'north_shift':north_shift, 'strike': strike, 'length': length} + # Print the beat.prior config file + + print(' priors:',file=file) + for prior in priors: + print(' %s: !beat.heart.Parameter' % prior,file=file) + print(' name: %s' % prior,file=file) + print(' form: Uniform',file=file) + + print(' lower:',file=file) + if prior in default_bounds: + for z in range(len(east_shift)): + print(' - %.2f' % default_bounds[prior][0],file=file) + elif prior in bounds: + for z in range(len(east_shift)): + print(' - %.2f' % bounds[prior][z],file=file) + else: + print('Error') + + print(' upper:',file=file) + if prior in default_bounds: + for z in range(len(east_shift)): + print(' - %.2f' % default_bounds[prior][1],file=file) + elif prior in bounds: + for z in range(len(east_shift)): + print(' - %.2f' % bounds[prior][z],file=file) + else: + print('Error') + + print(' testvalue:',file=file) + if prior in default_bounds: + for z in range(len(east_shift)): + print(' - %.2f' % default_bounds[prior][2],file=file) + elif prior in bounds: + for z in range(len(east_shift)): + print(' - %.2f' % bounds[prior][z],file=file) + else: + print('Error') + +######################################################################33 +def create_parser(): + parser = argparse.ArgumentParser(description=f'Get BEAT prior config for multiple fault segments from shapefile', + formatter_class=argparse.RawTextHelpFormatter) + + parser.add_argument('file', type=str, help='shapefile of fault segments, e.g. "fault.shp"') + parser.add_argument('-elon', '--event_lon', dest='event_lon', type=str, required=True, + help='event longitude') + parser.add_argument('-elat', '--event_lat', dest='event_lat', type=str, required=True, + help='event latitude') + return parser + +def cmd_line_parser(iargs=None): + parser = create_parser() + inps = parser.parse_args(args=iargs) + + return inps + +############################################################################## +def main(iargs=None): + inps = cmd_line_parser(iargs) + + # Get fault segements + event = [float(inps.event_lon), float(inps.event_lat)] + fault = getFaultSegments(event, inps.file) + + # Export BEAT config file + mfaultToBeat(fault) + print('Writing BEAT prior config in beatPriors_multiSegment.txt') + + +############################################################################## +if __name__ == '__main__': + main(sys.argv[1:]) diff --git a/beat/config.py b/beat/config.py index 1111f27a..5ef37465 100644 --- a/beat/config.py +++ b/beat/config.py @@ -298,6 +298,14 @@ class NonlinearGFConfig(GFConfig): default=True, help='Flag, for replacing the crust from the earthmodel' 'with crust from the crust2 model.') + use_epcrust = Bool.T( + default=True, + help='Flag, for replacing the crust from the earthmodel' + 'with crust from the EPcrust model.') + epcrust_radius = Float.T( + default=50, + help='Get the average EPcrust profile in the defined radius around Reference location' + 'radius is in km') replace_water = Bool.T( default=True, help='Flag, for replacing water layers in the crust2 model.') @@ -1818,6 +1826,8 @@ def init_config(name, date=None, min_magnitude=6.0, main_path='./', c.geodetic_config.gf_config.custom_velocity_model = \ load_model().extract(depth_max=100. * km) c.geodetic_config.gf_config.use_crust2 = False + c.geodetic_config.gf_config.use_epcrust = False + c.geodetic_config.gf_config.epcrust_radius = 0 c.geodetic_config.gf_config.replace_water = False else: c.geodetic_config = None @@ -1837,6 +1847,8 @@ def init_config(name, date=None, min_magnitude=6.0, main_path='./', c.seismic_config.gf_config.custom_velocity_model = \ load_model().extract(depth_max=100. * km) c.seismic_config.gf_config.use_crust2 = False + c.geodetic_config.gf_config.use_epcrust = False + c.geodetic_config.gf_config.epcrust_radius = 0 c.seismic_config.gf_config.replace_water = False else: c.seismic_config = None diff --git a/beat/heart.py b/beat/heart.py index eedfe603..696c788a 100644 --- a/beat/heart.py +++ b/beat/heart.py @@ -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 from pyrocko.cake import GradientLayer from pyrocko.fomosto import qseis, qssp from pyrocko.model import gnss @@ -1606,6 +1607,50 @@ def get_velocity_model( source_model = cake.load_model( earth_model_name, crust2_profile=profile) + +####################################################################################### + elif gfc.use_epcrust: + logger.info('Using EPcrust profile') + # load velocity profile from EPcrust0.5x0.5 and check for water + if gfc.epcrust_radius != 0: + profile = ep_crust.avg_profile(location.lat, location.lon, gfc.epcrust_radius) + logger.info('Average EPcrust profile in radius of %.f km around the coordinates: %.2f, %.2f' % (gfc.epcrust_radius, profile._lat, profile._lon)) + 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: + logger.debug('Replacing water layers! ...') + thickness_lwater = profile.get_layer(ep_crust.LWATER)[0] + + if thickness_lwater > 0.0: + logger.info( + 'Water layer %f in EPcrust model!' + ' Remove and add to lower crust' % thickness_lwater) + thickness_llowercrust = profile.get_layer( + ep_crust.LLOWERCRUST)[0] + thickness_lsed = profile.get_layer( + ep_crust.LSED)[0] + + profile.set_layer_thickness(ep_crust.LWATER, 0.0) + profile.set_layer_thickness( + ep_crust.LSED, + num.ceil(thickness_lsed / 3)) + profile.set_layer_thickness( + ep_crust.LLOWERCRUST, + thickness_llowercrust + + thickness_lwater + + (thickness_lsed - num.ceil(thickness_lsed / 3))) + + profile._elevation = 0.0 + logger.info('New Lower crust layer thickness %f' % + profile.get_layer(ep_crust.LLOWERCRUST)[0]) + else: + logger.debug('Not replacing water layers') + + source_model = cake.load_model( + earth_model_name, epcrust_profile=profile) +######################################################################################## + else: logger.info('Using global model ...') @@ -1713,16 +1758,16 @@ def get_fomosto_baseconfig( id='any_P', definition=definition)) - if 'any_S' in waveforms: + if 'any_S' in waveforms: # always add any_S to have automatic adjustment of GF window tabulated_phases.append(gf.TPDef( id='any_S', definition='s,S,s\\,S\\')) # surface waves - if 'slowest' in waveforms: - tabulated_phases.append(gf.TPDef( - id='slowest', - definition='0.8')) + #if 'slowest' in waveforms: + tabulated_phases.append(gf.TPDef( + id='slowest', + definition='0.8')) return gf.ConfigTypeA( id='%s_%s_%.3fHz_%s' % ( @@ -1840,20 +1885,38 @@ def choose_backend( pids = ['stored:' + wave for wave in waveforms] +################################################################################# + pids1 = pids[:-1] + pids2 = [pids[-1]] #slowest + ''' conf.time_region = ( gf.Timing( phase_defs=pids, offset=(-1.1 * window_extension), select='first'), gf.Timing( phase_defs=pids, offset=(1.6 * window_extension), select='last')) - + ''' + #MG + conf.time_region = ( + gf.Timing( + phase_defs=pids1, offset=(-1.2 * window_extension), select='first'), + gf.Timing( + phase_defs=pids2, offset=(2 * window_extension), select='last')) + ''' conf.cut = ( gf.Timing( phase_defs=pids, offset=-window_extension, select='first'), gf.Timing( phase_defs=pids, offset=(1.5 * window_extension), select='last')) + ''' + #MG + conf.cut = ( + gf.Timing( + phase_defs=pids1, offset=-window_extension, select='first'), + gf.Timing( + phase_defs=pids2, offset=(1.5 * window_extension), select='last')) conf.relevel_with_fade_in = True - + ''' conf.fade = ( gf.Timing( phase_defs=pids, offset=-window_extension, select='first'), @@ -1863,7 +1926,20 @@ def choose_backend( phase_defs=pids, offset=(window_extension), select='last'), gf.Timing( phase_defs=pids, offset=(1.6 * window_extension), select='last')) + ''' + #MG + conf.fade = ( + gf.Timing( + phase_defs=pids1, offset=-window_extension, select='first'), + gf.Timing( + phase_defs=pids1, offset=(-0.1) * window_extension, select='first'), + gf.Timing( + phase_defs=pids2, offset=(window_extension), select='last'), + gf.Timing( + phase_defs=pids2, offset=(2 * window_extension), select='last')) + +################################################################################### return conf