From dea687caca19e4bc13961f0f5ea271e19f7c82e7 Mon Sep 17 00:00:00 2001 From: James-Ravenell Date: Mon, 23 May 2022 10:22:59 -0600 Subject: [PATCH] formatted with black and flake 8 --- strec/cmt.py | 25 ++- strec/gcmt.py | 114 ++++++++---- strec/slab.py | 39 ++-- strec/subduction.py | 54 +++--- strec/subtype.py | 423 ++++++++++++++++++++++++-------------------- strec/tensor.py | 117 ++++++------ 6 files changed, 437 insertions(+), 335 deletions(-) diff --git a/strec/cmt.py b/strec/cmt.py index f3b7c86..581c247 100755 --- a/strec/cmt.py +++ b/strec/cmt.py @@ -14,7 +14,8 @@ def getComposite(rows): """Calculate composite moment tensor. Args: - rows (list):List of tuples containing (mrr,mtt,mpp,mrt,mrp,mtp) + rows (list): + List of tuples containing (mrr,mtt,mpp,mrt,mrp,mtp) Returns: tuple: Tuple of (composite moment tensor dictionary (see fill_tensor_from_angles), @@ -45,8 +46,7 @@ def getComposite(rows): vm13 = comp_sq_mean_sq[3] vm23 = comp_sq_mean_sq[4] vm12 = comp_sq_mean_sq[5] - varforbenius = vm11 * vm11 + vm12 * vm12 + \ - vm13 * vm13 + vm22 * vm22 + vm23 * vm23 + varforbenius = vm11 * vm11 + vm12 * vm12 + vm13 * vm13 + vm22 * vm22 + vm23 * vm23 forbenius = m11 * m11 + m12 * m12 + m13 * m13 + m22 * m22 + m23 * m23 similarity = np.sqrt(varforbenius) / forbenius @@ -61,8 +61,9 @@ def getComposite(rows): return (tensor, similarity, nrows) -def getCompositeCMT(lat, lon, depth, dbfile, box=0.1, depthbox=10, nmin=3, maxbox=1.0, - dbox=0.09): +def getCompositeCMT( + lat, lon, depth, dbfile, box=0.1, depthbox=10, nmin=3, maxbox=1.0, dbox=0.09 +): """Search a database for list of moment tensors, calculate composite moment tensor. Args: lat (float): Latitude (dd). @@ -85,10 +86,16 @@ def getCompositeCMT(lat, lon, depth, dbfile, box=0.1, depthbox=10, nmin=3, maxbo rows = [] searchwidth = box while len(rows) < nmin and searchwidth < maxbox: - qstr = ("SELECT mrr,mtt,mpp,mrt,mrp,mtp FROM earthquake WHERE lat >= %.4f AND " - "lat <= %.4f AND lon >= %.4f AND lon <= %.4f") - query = qstr % (lat - searchwidth, lat + searchwidth, - lon - searchwidth, lon + searchwidth) + qstr = ( + "SELECT mrr,mtt,mpp,mrt,mrp,mtp FROM earthquake WHERE lat >= %.4f AND " + "lat <= %.4f AND lon >= %.4f AND lon <= %.4f" + ) + query = qstr % ( + lat - searchwidth, + lat + searchwidth, + lon - searchwidth, + lon + searchwidth, + ) cursor.execute(query) rows = cursor.fetchall() if len(rows) >= nmin: diff --git a/strec/gcmt.py b/strec/gcmt.py index 046603f..9b198f8 100644 --- a/strec/gcmt.py +++ b/strec/gcmt.py @@ -12,8 +12,12 @@ TIMEOUT = 30 -HIST_GCMT_URL = 'http://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/jan76_dec17.ndk.gz' -MONTHLY_GCMT_URL = 'http://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/NEW_MONTHLY/' +HIST_GCMT_URL = ( + "http://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/jan76_dec17.ndk.gz" +) +MONTHLY_GCMT_URL = ( + "http://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/NEW_MONTHLY/" +) def fetch_gcmt(): @@ -34,13 +38,13 @@ def fetch_gcmt(): - mtp Mtp moment tensor component """ t1 = str(datetime.utcnow()) - print('%s - Fetching historical GCMT data...' % t1) + print("%s - Fetching historical GCMT data..." % t1) histfile = get_historical_gcmt() t2 = str(datetime.utcnow()) - print('%s - Converting to dataframe...' % t2) + print("%s - Converting to dataframe..." % t2) dataframe = ndk_to_dataframe(histfile) t3 = str(datetime.utcnow()) - print('%s - Fetching monthly data...' % t3) + print("%s - Fetching monthly data..." % t3) start_year = 2018 end_year = datetime.utcnow().year end_month = datetime.utcnow().month @@ -48,10 +52,10 @@ def fetch_gcmt(): for month in range(1, 13): if year == end_year and month > end_month: continue - print('Fetching GCMT data for %i month %i...' % (year, month)) + print("Fetching GCMT data for %i month %i..." % (year, month)) monthfile = get_monthly_gcmt(year, month) if monthfile is None: - print('No NDK file for %i month %i' % (year, month)) + print("No NDK file for %i month %i" % (year, month)) continue month_frame = ndk_to_dataframe(monthfile) dataframe = dataframe.append(month_frame) @@ -76,11 +80,11 @@ def get_historical_gcmt(): fh.close() handle, zipfile = tempfile.mkstemp() os.close(handle) - f = open(zipfile, 'wb') + f = open(zipfile, "wb") f.write(zip_bytes) f.close() - gz = gzip.GzipFile(zipfile, mode='rb') - unzip_bytes = gz.read().decode('utf-8') + gz = gzip.GzipFile(zipfile, mode="rb") + unzip_bytes = gz.read().decode("utf-8") string_unzip = io.StringIO(unzip_bytes) return string_unzip except Exception as msg: @@ -102,14 +106,27 @@ def get_monthly_gcmt(year, month): io.StringIO: StringIO object containing NDK file for given month/year. """ - strmonth = ['jan', 'feb', 'mar', 'apr', 'may', 'jun', - 'jul', 'aug', 'sep', 'oct', 'nov', 'dec'][month - 1] + strmonth = [ + "jan", + "feb", + "mar", + "apr", + "may", + "jun", + "jul", + "aug", + "sep", + "oct", + "nov", + "dec", + ][month - 1] stryear = str(year) - url = parse.urljoin(MONTHLY_GCMT_URL, os.path.join( - stryear, strmonth + stryear[2:] + '.ndk')) + url = parse.urljoin( + MONTHLY_GCMT_URL, os.path.join(stryear, strmonth + stryear[2:] + ".ndk") + ) try: fh = request.urlopen(url, timeout=TIMEOUT) - bytes = fh.read().decode('utf-8') + bytes = fh.read().decode("utf-8") fh.close() stringfile = io.StringIO(bytes) return stringfile @@ -138,12 +155,25 @@ def ndk_to_dataframe(ndkfile): - mrp Mrp moment tensor component - mtp Mtp moment tensor component """ - if not hasattr(ndkfile, 'read'): - ndkfile = open(ndkfile, 'r') + if not hasattr(ndkfile, "read"): + ndkfile = open(ndkfile, "r") lc = 0 - df = pd.DataFrame(columns=['time', 'lat', 'lon', 'depth', 'mag', - 'mrr', 'mtt', 'mpp', 'mrt', 'mrp', 'mtp']) + df = pd.DataFrame( + columns=[ + "time", + "lat", + "lon", + "depth", + "mag", + "mrr", + "mtt", + "mpp", + "mrt", + "mrp", + "mtp", + ] + ) tdict = {} for line in ndkfile.readlines(): if (lc + 1) % 5 == 1: @@ -163,7 +193,7 @@ def ndk_to_dataframe(ndkfile): if (lc + 1) % 5 == 0: _parse_line5(line, tdict) lc += 1 - tdict.pop('exponent') # remove now extraneous field + tdict.pop("exponent") # remove now extraneous field df = df.append(tdict, ignore_index=True) tdict = {} continue @@ -175,6 +205,12 @@ def ndk_to_dataframe(ndkfile): def _parse_line1(line, tdict): """Parse the first line of an NDK file. + Args: + line (str/byte): + The line of the NDK file + tdict (dict): + dictionary to inpiut values from ndk files + """ dstr = line[5:26] year = int(dstr[0:4]) @@ -190,31 +226,39 @@ def _parse_line1(line, tdict): if microseconds > 999999: microseconds = 999999 - tdict['time'] = datetime(year, month, day, hour, - minute, seconds, microseconds) + tdict["time"] = datetime(year, month, day, hour, minute, seconds, microseconds) - tdict['lat'] = float(line[27:33]) - tdict['lon'] = float(line[34:41]) - tdict['depth'] = float(line[42:47]) + tdict["lat"] = float(line[27:33]) + tdict["lon"] = float(line[34:41]) + tdict["depth"] = float(line[42:47]) def _parse_line4(line, tdict): """Parse the fourth line of an NDK file. + Args: + line (str/byte): + The line of the NDK file + tdict (dict): + dictionary to inpiut values from ndk files """ - tdict['exponent'] = float(line[0:2]) - tdict['mrr'] = float(line[2:9]) * np.power(10.0, tdict['exponent']) - tdict['mtt'] = float(line[15:22]) * np.power(10.0, tdict['exponent']) - tdict['mpp'] = float(line[28:35]) * np.power(10.0, tdict['exponent']) - tdict['mrt'] = float(line[41:48]) * np.power(10.0, tdict['exponent']) - tdict['mrp'] = float(line[54:61]) * np.power(10.0, tdict['exponent']) - tdict['mtp'] = float(line[67:74]) * np.power(10.0, tdict['exponent']) + tdict["exponent"] = float(line[0:2]) + tdict["mrr"] = float(line[2:9]) * np.power(10.0, tdict["exponent"]) + tdict["mtt"] = float(line[15:22]) * np.power(10.0, tdict["exponent"]) + tdict["mpp"] = float(line[28:35]) * np.power(10.0, tdict["exponent"]) + tdict["mrt"] = float(line[41:48]) * np.power(10.0, tdict["exponent"]) + tdict["mrp"] = float(line[54:61]) * np.power(10.0, tdict["exponent"]) + tdict["mtp"] = float(line[67:74]) * np.power(10.0, tdict["exponent"]) def _parse_line5(line, tdict): """Parse the fifth line of an NDK file. + Args: + line (str/byte): + The line of the NDK file + tdict (dict): + dictionary to inpiut values from ndk files """ - scalar_moment = float(line[49:56].strip()) * \ - np.power(10.0, tdict['exponent']) - tdict['mag'] = ((2.0 / 3.0) * np.log10(scalar_moment)) - 10.7 + scalar_moment = float(line[49:56].strip()) * np.power(10.0, tdict["exponent"]) + tdict["mag"] = ((2.0 / 3.0) * np.log10(scalar_moment)) - 10.7 diff --git a/strec/slab.py b/strec/slab.py index fb01b98..83e3e7b 100644 --- a/strec/slab.py +++ b/strec/slab.py @@ -14,8 +14,7 @@ class GridSlab(object): - """Represents USGS Slab model grids for a given subduction zone. - """ + """Represents USGS Slab model grids for a given subduction zone.""" def __init__(self, depth_file, dip_file, strike_file, error_file): """Construct GridSlab object from input files. @@ -35,7 +34,7 @@ def __init__(self, depth_file, dip_file, strike_file, error_file): # as all of the slab grids. Read it into a local dictionary if found, # otherwise we'll use the MAX_INTERFACE_DEPTH constant found above. fpath, fname = os.path.split(self._depth_file) - table_file_name = os.path.join(fpath, 'maximum_interface_depths.csv') + table_file_name = os.path.join(fpath, "maximum_interface_depths.csv") if os.path.isfile(table_file_name): self._slab_table = pd.read_csv(table_file_name) else: @@ -70,7 +69,7 @@ def getSlabInfo(self, lat, lon): lat (float): Hypocentral latitude in decimal degrees. lon (float): Hypocentral longitude in decimal degrees. Returns: - dict: Dictionary containing keys: + slabinfo (dict): Dictionary containing keys: - region Three letter Slab model region code. - strike Slab model strike angle. - dip Slab model dip angle. @@ -81,7 +80,7 @@ def getSlabInfo(self, lat, lon): if not self.contains(lat, lon): return slabinfo fpath, fname = os.path.split(self._depth_file) - parts = fname.split('_') + parts = fname.split("_") region = parts[0] depth_grid = GMTGrid.load(self._depth_file) # slab grids are negative depth @@ -109,16 +108,18 @@ def getSlabInfo(self, lat, lon): # get the maximum interface depth from table (if present) if self._slab_table is not None: df = self._slab_table - max_int_depth = df[df['zone'] == region].iloc[0]['interface_max_depth'] + max_int_depth = df[df["zone"] == region].iloc[0]["interface_max_depth"] else: max_int_depth = MAX_INTERFACE_DEPTH - slabinfo = {'region': region, - 'strike': strike, - 'dip': dip, - 'depth': depth, - 'maximum_interface_depth' : max_int_depth, - 'depth_uncertainty': error} + slabinfo = { + "region": region, + "strike": strike, + "dip": dip, + "depth": depth, + "maximum_interface_depth": max_int_depth, + "depth_uncertainty": error, + } return slabinfo @@ -132,7 +133,7 @@ def __init__(self, datafolder): Args: datafolder (str): String path where grid files and GeoJSON file reside. """ - self._depth_files = glob.glob(os.path.join(datafolder, '*_dep*.grd')) + self._depth_files = glob.glob(os.path.join(datafolder, "*_dep*.grd")) def getSlabInfo(self, lat, lon, depth): """Query the entire set of slab models and return a SlabInfo object, or None. @@ -143,7 +144,7 @@ def getSlabInfo(self, lat, lon, depth): depth (float): Hypocentral depth in km. Returns: - dict: Dictionary containing keys: + slabinfo (dict): Dictionary containing keys: - region Three letter Slab model region code. - strike Slab model strike angle. - dip Slab model dip angle. @@ -155,9 +156,9 @@ def getSlabInfo(self, lat, lon, depth): slabinfo = {} # loop over all slab regions, return keep all slabs found for depth_file in self._depth_files: - dip_file = depth_file.replace('dep', 'dip') - strike_file = depth_file.replace('dep', 'str') - error_file = depth_file.replace('dep', 'unc') + dip_file = depth_file.replace("dep", "dip") + strike_file = depth_file.replace("dep", "str") + error_file = depth_file.replace("dep", "unc") if not os.path.isfile(error_file): error_file = None gslab = GridSlab(depth_file, dip_file, strike_file, error_file) @@ -165,11 +166,11 @@ def getSlabInfo(self, lat, lon, depth): if not len(tslabinfo): continue else: - depth = tslabinfo['depth'] + depth = tslabinfo["depth"] if depth < deep_depth: slabinfo = tslabinfo.copy() deep_depth = depth - elif np.isnan(depth) and 'depth' not in slabinfo: + elif np.isnan(depth) and "depth" not in slabinfo: slabinfo = tslabinfo.copy() return slabinfo diff --git a/strec/subduction.py b/strec/subduction.py index 0c9e3c8..e4dcf05 100644 --- a/strec/subduction.py +++ b/strec/subduction.py @@ -27,11 +27,11 @@ def __init__(self, slab_params, tensor_params, depth, config): event. DDEPTH_INTRA - intra-slab depth range. """ - self._dstrike = float(config['CONSTANTS']['dstrike_interf']) - self._ddip = float(config['CONSTANTS']['ddip_interf']) - self._dlambda = float(config['CONSTANTS']['dlambda']) - self._ddepth_interface = float(config['CONSTANTS']['ddepth_interf']) - self._ddepth_intraslab = float(config['CONSTANTS']['ddepth_intra']) + self._dstrike = float(config["CONSTANTS"]["dstrike_interf"]) + self._ddip = float(config["CONSTANTS"]["ddip_interf"]) + self._dlambda = float(config["CONSTANTS"]["dlambda"]) + self._ddepth_interface = float(config["CONSTANTS"]["ddepth_interf"]) + self._ddepth_intraslab = float(config["CONSTANTS"]["ddepth_intra"]) self._slab_params = slab_params.copy() if isinstance(tensor_params, dict): self._tensor_params = tensor_params.copy() @@ -42,18 +42,17 @@ def __init__(self, slab_params, tensor_params, depth, config): def checkRupturePlane(self): """Compare moment tensor angles to slab angles, return True if similar. - Returns: bool: Boolean value indicating if moment tensor is similar to slab. """ if self._tensor_params is None: return False - strike = self._slab_params['strike'] - 90 - a = self._tensor_params['P']['azimuth'] - b1 = (norm_angle(strike) - self._dstrike) - b2 = (norm_angle(strike) + self._dstrike) - b3 = (norm_angle(strike) - self._dstrike) - b4 = (norm_angle(strike) + self._dstrike) + strike = self._slab_params["strike"] - 90 + a = self._tensor_params["P"]["azimuth"] + b1 = norm_angle(strike) - self._dstrike + b2 = norm_angle(strike) + self._dstrike + b3 = norm_angle(strike) - self._dstrike + b4 = norm_angle(strike) + self._dstrike if a > 270 and b1 < 90: b1 = b1 + 360 @@ -79,15 +78,18 @@ def checkRupturePlane(self): c4 = a <= b4 m2a = (c1 and c2) or (c3 and c4) - m2b = ((self._tensor_params['P']['plunge'] >= self._slab_params['dip'] - - self._ddip) and - (self._tensor_params['P']['plunge'] <= self._slab_params['dip'] + - self._ddip)) - m2c1 = ((self._tensor_params['NP1']['rake'] > 90 - self._dlambda) and - (self._tensor_params['NP1']['rake'] < 90 + self._dlambda)) - m2c2 = ((self._tensor_params['NP2']['rake'] > 90 - self._dlambda) and - (self._tensor_params['NP2']['rake'] < 90 + self._dlambda)) - if (m2a and m2b and (m2c1 or m2c2)): + m2b = ( + self._tensor_params["P"]["plunge"] >= self._slab_params["dip"] - self._ddip + ) and ( + self._tensor_params["P"]["plunge"] <= self._slab_params["dip"] + self._ddip + ) + m2c1 = (self._tensor_params["NP1"]["rake"] > 90 - self._dlambda) and ( + self._tensor_params["NP1"]["rake"] < 90 + self._dlambda + ) + m2c2 = (self._tensor_params["NP2"]["rake"] > 90 - self._dlambda) and ( + self._tensor_params["NP2"]["rake"] < 90 + self._dlambda + ) + if m2a and m2b and (m2c1 or m2c2): return True else: return False @@ -98,9 +100,9 @@ def checkInterfaceDepth(self): Returns: bool: True if event is close to slab interface depth, False otherwise. """ - c1 = self._depth >= self._slab_params['depth'] - self._ddepth_interface - c2 = self._depth < self._slab_params['depth'] + self._ddepth_interface - if (c1 and c2): + c1 = self._depth >= self._slab_params["depth"] - self._ddepth_interface + c2 = self._depth < self._slab_params["depth"] + self._ddepth_interface + if c1 and c2: return True else: return False @@ -111,9 +113,9 @@ def checkSlabDepth(self, intraslabdepth): Args: intraslabdepth (float): Upper limit of intraslab events. - @return: True if depth is deeper than the slab, False if not. + return: True if depth is deeper than the slab, False if not. """ - c1 = self._depth >= self._slab_params['depth'] - self._ddepth_intraslab + c1 = self._depth >= self._slab_params["depth"] - self._ddepth_intraslab c2 = self._depth >= intraslabdepth if c1 and c2: return True diff --git a/strec/subtype.py b/strec/subtype.py index 1aee8ea..7c3e6b4 100644 --- a/strec/subtype.py +++ b/strec/subtype.py @@ -16,55 +16,58 @@ from strec.tensor import fill_tensor_from_components from strec.utils import get_config -EVENT_URL = 'https://earthquake.usgs.gov/fdsnws/event/1/query?eventid=EVENTID&format=geojson' +EVENT_URL = ( + "https://earthquake.usgs.gov/fdsnws/event/1/query?eventid=EVENTID&format=geojson" +) SLAB_RAKE = 90 # presumed rake angle of slabs -SLAB_REGIONS = {'alu': 'Alaska-Aleutians', - 'cal': 'Calabria', - 'cam': 'Central America', - 'car': 'Caribbean', - 'cas': 'Cascadia', - 'cot': 'Cotabato', - 'hal': 'Halmahera', - 'hel': 'Helanic', - 'him': 'Himalaya', - 'hin': 'Hindu Kush', - 'izu': 'Izu-Bonin', - 'ker': 'Kermadec-Tonga', - 'kur': 'Kamchatka/Kurils/Japan', - 'mak': 'Makran', - 'man': 'Manila', - 'mex': 'Central America', - 'mue': 'Muertos', - 'pam': 'Pamir', - 'pan': 'Panama', - 'phi': 'Philippines', - 'png': 'New Guinea', - 'puy': 'Puysegur', - 'ryu': 'Ryukyu', - 'sam': 'South America', - 'sco': 'Scotia', - 'sol': 'Solomon Islands', - 'sul': 'Sulawesi', - 'sum': 'Sumatra-Java', - 'van': 'Santa Cruz Islands/Vanuatu/Loyalty Islands'} - -CONSTANTS = {'tplunge_rs': 50, - 'bplunge_ds': 30, - 'bplunge_ss': 55, - 'pplunge_nm': 55, - 'delplunge_ss': 20} +SLAB_REGIONS = { + "alu": "Alaska-Aleutians", + "cal": "Calabria", + "cam": "Central America", + "car": "Caribbean", + "cas": "Cascadia", + "cot": "Cotabato", + "hal": "Halmahera", + "hel": "Helanic", + "him": "Himalaya", + "hin": "Hindu Kush", + "izu": "Izu-Bonin", + "ker": "Kermadec-Tonga", + "kur": "Kamchatka/Kurils/Japan", + "mak": "Makran", + "man": "Manila", + "mex": "Central America", + "mue": "Muertos", + "pam": "Pamir", + "pan": "Panama", + "phi": "Philippines", + "png": "New Guinea", + "puy": "Puysegur", + "ryu": "Ryukyu", + "sam": "South America", + "sco": "Scotia", + "sol": "Solomon Islands", + "sul": "Sulawesi", + "sum": "Sumatra-Java", + "van": "Santa Cruz Islands/Vanuatu/Loyalty Islands", +} + +CONSTANTS = { + "tplunge_rs": 50, + "bplunge_ds": 30, + "bplunge_ss": 55, + "pplunge_nm": 55, + "delplunge_ss": 20, +} class SubductionSelector(object): - """For events that are inside a subduction zone, determine subduction zone properties. - - """ + """For events that are inside a subduction zone, determine subduction zone + properties.""" def __init__(self, prefix=None, verbose=False): - """Construct a SubductionSelector object. - - """ + """Construct a SubductionSelector object.""" if prefix is not None: self.logger = logging.getLogger(prefix) else: @@ -79,7 +82,8 @@ def getSubductionTypeByID(self, eventid): Args: eventid (str): ComCat EventID (Sumatra is official20041226005853450_30). Returns: - Pandas Series object with indices: + + results (Panda series): Pandas Series object with indices: - TectonicRegion : (Subduction,Active,Stable,Volcanic) - TectonicDomain : SZ (generic) - FocalMechanism : (RS [Reverse],SS [Strike-Slip], NM [Normal], ALL @@ -120,17 +124,16 @@ def getSubductionTypeByID(self, eventid): AttributeError if the eventid is not found in ComCat. """ if self.verbose: - self.logger.info('Inside getSubductionTypeByID...') + self.logger.info("Inside getSubductionTypeByID...") lat, lon, depth, tensor_params = self.getOnlineTensor(eventid) if self.verbose: - self.logger.info('Tensor Parameters: %s' % str(tensor_params)) + self.logger.info("Tensor Parameters: %s" % str(tensor_params)) if lat is None: - raise AttributeError('Event %s is not found in ComCat.' % eventid) + raise AttributeError("Event %s is not found in ComCat." % eventid) lat = float(lat) lon = float(lon) - results = self.getSubductionType( - lat, lon, depth, tensor_params=tensor_params) + results = self.getSubductionType(lat, lon, depth, tensor_params=tensor_params) return results def getOnlineTensor(self, eventid): @@ -162,7 +165,7 @@ def getOnlineTensor(self, eventid): - rake """ if self.verbose: - self.logger.info('Inside getOnlineTensor') + self.logger.info("Inside getOnlineTensor") try: detail = get_event_by_id(eventid) except Exception as e: @@ -173,97 +176,102 @@ def getOnlineTensor(self, eventid): lat = detail.latitude lon = detail.longitude depth = detail.depth - if not detail.hasProduct('moment-tensor'): - self.logger.info('No moment tensor available for %s' % eventid) + if not detail.hasProduct("moment-tensor"): + self.logger.info("No moment tensor available for %s" % eventid) return lat, lon, depth, None if self.verbose: - self.logger.info('Getting tensor components...') - tensor = detail.getProducts('moment-tensor')[0] + self.logger.info("Getting tensor components...") + tensor = detail.getProducts("moment-tensor")[0] tensor_params = {} - btype = 'unknown' - if tensor.hasProperty('derived-magnitude-type'): - btype = tensor['derived-magnitude-type'] - elif tensor.hasProperty('beachball-type'): - btype = tensor['beachball-type'] - if btype.find('/') > -1: - btype = btype.split('/')[-1] - tensor_params['type'] = btype - tensor_params['source'] = tensor['eventsource'] + \ - '_' + tensor['eventsourcecode'] - - tensor_params['mtt'] = float(tensor['tensor-mtt']) - tensor_params['mpp'] = float(tensor['tensor-mpp']) - tensor_params['mrr'] = float(tensor['tensor-mrr']) - tensor_params['mtp'] = float(tensor['tensor-mtp']) - tensor_params['mrt'] = float(tensor['tensor-mrt']) - tensor_params['mrp'] = float(tensor['tensor-mrp']) + btype = "unknown" + if tensor.hasProperty("derived-magnitude-type"): + btype = tensor["derived-magnitude-type"] + elif tensor.hasProperty("beachball-type"): + btype = tensor["beachball-type"] + if btype.find("/") > -1: + btype = btype.split("/")[-1] + tensor_params["type"] = btype + tensor_params["source"] = ( + tensor["eventsource"] + "_" + tensor["eventsourcecode"] + ) + + tensor_params["mtt"] = float(tensor["tensor-mtt"]) + tensor_params["mpp"] = float(tensor["tensor-mpp"]) + tensor_params["mrr"] = float(tensor["tensor-mrr"]) + tensor_params["mtp"] = float(tensor["tensor-mtp"]) + tensor_params["mrt"] = float(tensor["tensor-mrt"]) + tensor_params["mrp"] = float(tensor["tensor-mrp"]) if self.verbose: - self.logger.info('Getting tensor axes...') + self.logger.info("Getting tensor axes...") # sometimes the online MT is missing properties - if not tensor.hasProperty('t-axis-length'): + if not tensor.hasProperty("t-axis-length"): if self.verbose: - self.logger.info('Calling fill_tensor function...') - tensor_dict = fill_tensor_from_components(tensor_params['mrr'], - tensor_params['mtt'], - tensor_params['mpp'], - tensor_params['mrt'], - tensor_params['mrp'], - tensor_params['mtp']) - tensor_params['T'] = tensor_dict['T'].copy() - tensor_params['N'] = tensor_dict['T'].copy() - tensor_params['P'] = tensor_dict['P'].copy() + self.logger.info("Calling fill_tensor function...") + tensor_dict = fill_tensor_from_components( + tensor_params["mrr"], + tensor_params["mtt"], + tensor_params["mpp"], + tensor_params["mrt"], + tensor_params["mrp"], + tensor_params["mtp"], + ) + tensor_params["T"] = tensor_dict["T"].copy() + tensor_params["N"] = tensor_dict["T"].copy() + tensor_params["P"] = tensor_dict["P"].copy() else: T = {} - T['value'] = float(tensor['t-axis-length']) - T['plunge'] = float(tensor['t-axis-plunge']) - T['azimuth'] = float(tensor['t-axis-azimuth']) - tensor_params['T'] = T.copy() + T["value"] = float(tensor["t-axis-length"]) + T["plunge"] = float(tensor["t-axis-plunge"]) + T["azimuth"] = float(tensor["t-axis-azimuth"]) + tensor_params["T"] = T.copy() N = {} - N['value'] = float(tensor['n-axis-length']) - N['plunge'] = float(tensor['n-axis-plunge']) - N['azimuth'] = float(tensor['n-axis-azimuth']) - tensor_params['N'] = N.copy() + N["value"] = float(tensor["n-axis-length"]) + N["plunge"] = float(tensor["n-axis-plunge"]) + N["azimuth"] = float(tensor["n-axis-azimuth"]) + tensor_params["N"] = N.copy() P = {} - P['value'] = float(tensor['p-axis-length']) - P['plunge'] = float(tensor['p-axis-plunge']) - P['azimuth'] = float(tensor['p-axis-azimuth']) - tensor_params['P'] = P.copy() + P["value"] = float(tensor["p-axis-length"]) + P["plunge"] = float(tensor["p-axis-plunge"]) + P["azimuth"] = float(tensor["p-axis-azimuth"]) + tensor_params["P"] = P.copy() if self.verbose: - self.logger.info('Getting tensor angles...') - if not tensor.hasProperty('nodal-plane-1-strike'): + self.logger.info("Getting tensor angles...") + if not tensor.hasProperty("nodal-plane-1-strike"): if self.verbose: - self.logger.info('Calling fill_tensor function...') - tensor2 = fill_tensor_from_components(tensor_params['mrr'], - tensor_params['mtt'], - tensor_params['mpp'], - tensor_params['mrt'], - tensor_params['mrp'], - tensor_params['mtp']) - tensor_params['NP1'] = tensor2['NP1'].copy() - tensor_params['NP2'] = tensor2['NP2'].copy() + self.logger.info("Calling fill_tensor function...") + tensor2 = fill_tensor_from_components( + tensor_params["mrr"], + tensor_params["mtt"], + tensor_params["mpp"], + tensor_params["mrt"], + tensor_params["mrp"], + tensor_params["mtp"], + ) + tensor_params["NP1"] = tensor2["NP1"].copy() + tensor_params["NP2"] = tensor2["NP2"].copy() else: NP1 = {} - NP1['strike'] = float(tensor['nodal-plane-1-strike']) - NP1['dip'] = float(tensor['nodal-plane-1-dip']) - if 'nodal-plane-1-rake' in tensor.properties: - NP1['rake'] = float(tensor['nodal-plane-1-rake']) + NP1["strike"] = float(tensor["nodal-plane-1-strike"]) + NP1["dip"] = float(tensor["nodal-plane-1-dip"]) + if "nodal-plane-1-rake" in tensor.properties: + NP1["rake"] = float(tensor["nodal-plane-1-rake"]) else: - NP1['rake'] = float(tensor['nodal-plane-1-slip']) - tensor_params['NP1'] = NP1.copy() + NP1["rake"] = float(tensor["nodal-plane-1-slip"]) + tensor_params["NP1"] = NP1.copy() NP2 = {} - NP2['strike'] = float(tensor['nodal-plane-2-strike']) - NP2['dip'] = float(tensor['nodal-plane-2-dip']) - if 'nodal-plane-2-rake' in tensor.properties: - NP2['rake'] = float(tensor['nodal-plane-2-rake']) + NP2["strike"] = float(tensor["nodal-plane-2-strike"]) + NP2["dip"] = float(tensor["nodal-plane-2-dip"]) + if "nodal-plane-2-rake" in tensor.properties: + NP2["rake"] = float(tensor["nodal-plane-2-rake"]) else: - NP2['rake'] = float(tensor['nodal-plane-2-slip']) - tensor_params['NP2'] = NP2.copy() + NP2["rake"] = float(tensor["nodal-plane-2-slip"]) + tensor_params["NP2"] = NP2.copy() return lat, lon, depth, tensor_params @@ -298,7 +306,7 @@ def getSubductionType(self, lat, lon, depth, eventid=None, tensor_params=None): (optional) - source Moment Tensor source (regional network, name of study, etc.) Returns: - Pandas Series object with indices: + results (series):Pandas Series object with indices: - TectonicRegion : (Subduction,Active,Stable,Volcanic) - FocalMechanism : (RS [Reverse],SS [Strike-Slip], NM [Normal], ALL [Unknown]) @@ -325,14 +333,14 @@ def getSubductionType(self, lat, lon, depth, eventid=None, tensor_params=None): - SlabModelStrike : Strike of slab at epicenter. """ if self.verbose: - self.logger.info('Inside getSubductionType...') + self.logger.info("Inside getSubductionType...") # sometimes events are specified with negative depths, which don't work with # our algorithms. Pin those depths to 0. if depth < 0: depth = 0 config = self._config - slab_data_folder = config['DATA']['slabfolder'] + slab_data_folder = config["DATA"]["slabfolder"] tensor_type = None tensor_source = None similarity = np.nan @@ -341,90 +349,115 @@ def getSubductionType(self, lat, lon, depth, eventid=None, tensor_params=None): if eventid is not None: tlat, tlon, tdep, tensor_params = self.getOnlineTensor(eventid) if tensor_params is not None: - tensor_type = tensor_params['type'] - tensor_source = tensor_params['source'] + tensor_type = tensor_params["type"] + tensor_source = tensor_params["source"] if tensor_params is None: dbfile = os.path.join( - config['DATA']['folder'], config['DATA']['dbfile']) - minboxcomp = float(config['CONSTANTS']['minradial_distcomp']) - maxboxcomp = float(config['CONSTANTS']['maxradial_distcomp']) - dboxcomp = float(config['CONSTANTS']['step_distcomp']) - depthboxcomp = float(config['CONSTANTS']['depth_rangecomp']) + config["DATA"]["folder"], config["DATA"]["dbfile"] + ) + minboxcomp = float(config["CONSTANTS"]["minradial_distcomp"]) + maxboxcomp = float(config["CONSTANTS"]["maxradial_distcomp"]) + dboxcomp = float(config["CONSTANTS"]["step_distcomp"]) + depthboxcomp = float(config["CONSTANTS"]["depth_rangecomp"]) # Minimum number of events required to get composite mechanism - nmin = int(config['CONSTANTS']['minno_comp']) - tensor_params, similarity, nevents = getCompositeCMT(lat, lon, - depth, dbfile, - box=minboxcomp, - depthbox=depthboxcomp, - maxbox=maxboxcomp, - dbox=dboxcomp, - nmin=nmin) + nmin = int(config["CONSTANTS"]["minno_comp"]) + tensor_params, similarity, nevents = getCompositeCMT( + lat, + lon, + depth, + dbfile, + box=minboxcomp, + depthbox=depthboxcomp, + maxbox=maxboxcomp, + dbox=dboxcomp, + nmin=nmin, + ) if tensor_params is not None: - tensor_type = 'composite' - tensor_source = 'composite' + tensor_type = "composite" + tensor_source = "composite" else: - if 'type' in tensor_params: - tensor_type = tensor_params['type'] - if 'source' in tensor_params: - tensor_source = tensor_params['source'] + if "type" in tensor_params: + tensor_type = tensor_params["type"] + if "source" in tensor_params: + tensor_source = tensor_params["source"] slab_collection = SlabCollection(slab_data_folder) slab_params = slab_collection.getSlabInfo(lat, lon, depth) results = self._regionalizer.getRegions(lat, lon, depth) - results['TensorType'] = tensor_type - results['TensorSource'] = tensor_source - results['CompositeVariability'] = similarity - results['NComposite'] = nevents - results['FocalMechanism'] = get_focal_mechanism(tensor_params) + results["TensorType"] = tensor_type + results["TensorSource"] = tensor_source + results["CompositeVariability"] = similarity + results["NComposite"] = nevents + results["FocalMechanism"] = get_focal_mechanism(tensor_params) if len(slab_params): - if np.isnan(slab_params['depth']): - results['SlabModelRegion'] = SLAB_REGIONS[slab_params['region']] - results['KaganAngle'] = np.nan + if np.isnan(slab_params["depth"]): + results["SlabModelRegion"] = SLAB_REGIONS[slab_params["region"]] + results["KaganAngle"] = np.nan else: - results['SlabModelRegion'] = SLAB_REGIONS[slab_params['region']] - results['SlabModelDepth'] = slab_params['depth'] - results['SlabModelDepthUncertainty'] = slab_params['depth_uncertainty'] - results['SlabModelDip'] = slab_params['dip'] - results['SlabModelStrike'] = slab_params['strike'] - results['SlabModelMaximumDepth'] = slab_params['maximum_interface_depth' - ] + results["SlabModelRegion"] = SLAB_REGIONS[slab_params["region"]] + results["SlabModelDepth"] = slab_params["depth"] + results["SlabModelDepthUncertainty"] = slab_params["depth_uncertainty"] + results["SlabModelDip"] = slab_params["dip"] + results["SlabModelStrike"] = slab_params["strike"] + results["SlabModelMaximumDepth"] = slab_params[ + "maximum_interface_depth" + ] if tensor_params is not None: - np1 = tensor_params['NP1'] - kagan = get_kagan_angle(slab_params['strike'], slab_params['dip'], - SLAB_RAKE, np1['strike'], np1['dip'], - np1['rake']) - results['KaganAngle'] = kagan + np1 = tensor_params["NP1"] + kagan = get_kagan_angle( + slab_params["strike"], + slab_params["dip"], + SLAB_RAKE, + np1["strike"], + np1["dip"], + np1["rake"], + ) + results["KaganAngle"] = kagan else: - results['KaganAngle'] = np.nan + results["KaganAngle"] = np.nan else: - results['FocalMechanism'] = get_focal_mechanism(tensor_params) - results['SlabModelRegion'] = '' - results['SlabModelDepth'] = np.nan - results['SlabModelDepthUncertainty'] = np.nan - results['SlabModelDip'] = np.nan - results['SlabModelStrike'] = np.nan - results['SlabModelMaximumDepth'] = np.nan - results['KaganAngle'] = np.nan - - results = results.reindex(index=['TectonicRegion', 'FocalMechanism', - 'TensorType', 'TensorSource', 'KaganAngle', - 'CompositeVariability', 'NComposite', - 'DistanceToStable', 'DistanceToActive', - 'DistanceToSubduction', 'DistanceToVolcanic', - 'Oceanic', 'DistanceToOceanic', - 'DistanceToContinental', 'SlabModelRegion', - 'SlabModelDepth', 'SlabModelDepthUncertainty', - 'SlabModelDip', 'SlabModelStrike', - 'SlabModelMaximumDepth']) + results["FocalMechanism"] = get_focal_mechanism(tensor_params) + results["SlabModelRegion"] = "" + results["SlabModelDepth"] = np.nan + results["SlabModelDepthUncertainty"] = np.nan + results["SlabModelDip"] = np.nan + results["SlabModelStrike"] = np.nan + results["SlabModelMaximumDepth"] = np.nan + results["KaganAngle"] = np.nan + + results = results.reindex( + index=[ + "TectonicRegion", + "FocalMechanism", + "TensorType", + "TensorSource", + "KaganAngle", + "CompositeVariability", + "NComposite", + "DistanceToStable", + "DistanceToActive", + "DistanceToSubduction", + "DistanceToVolcanic", + "Oceanic", + "DistanceToOceanic", + "DistanceToContinental", + "SlabModelRegion", + "SlabModelDepth", + "SlabModelDepthUncertainty", + "SlabModelDip", + "SlabModelStrike", + "SlabModelMaximumDepth", + ] + ) return results def get_focal_mechanism(tensor_params): - """ Return focal mechanism (strike-slip,normal, or reverse). + """Return focal mechanism (strike-slip,normal, or reverse). Args: tensor_params (dict): Dictionary containing the following fields: @@ -446,20 +479,20 @@ def get_focal_mechanism(tensor_params): str: Focal mechanism string 'SS','RS','NM',or 'ALL'. """ if tensor_params is None: - return 'ALL' + return "ALL" # implement eq 1 here - Tp = tensor_params['T']['plunge'] - Np = tensor_params['N']['plunge'] - Pp = tensor_params['P']['plunge'] - tplunge_rs = CONSTANTS['tplunge_rs'] - bplunge_ds = CONSTANTS['bplunge_ds'] - bplunge_ss = CONSTANTS['bplunge_ss'] - pplunge_nm = CONSTANTS['pplunge_nm'] - delplunge_ss = CONSTANTS['delplunge_ss'] + Tp = tensor_params["T"]["plunge"] + Np = tensor_params["N"]["plunge"] + Pp = tensor_params["P"]["plunge"] + tplunge_rs = CONSTANTS["tplunge_rs"] + bplunge_ds = CONSTANTS["bplunge_ds"] + bplunge_ss = CONSTANTS["bplunge_ss"] + pplunge_nm = CONSTANTS["pplunge_nm"] + delplunge_ss = CONSTANTS["delplunge_ss"] if Tp >= tplunge_rs and Np <= bplunge_ds: - return 'RS' + return "RS" if Np >= bplunge_ss and (Tp >= Pp - delplunge_ss and Tp <= Pp + delplunge_ss): - return 'SS' + return "SS" if Pp >= pplunge_nm and Np <= bplunge_ds: - return 'NM' - return 'ALL' + return "NM" + return "ALL" diff --git a/strec/tensor.py b/strec/tensor.py index 0b95b55..30ca63e 100644 --- a/strec/tensor.py +++ b/strec/tensor.py @@ -3,8 +3,9 @@ from obspy.imaging.beachball import aux_plane, mt2axes, mt2plane, MomentTensor -def fill_tensor_from_angles(strike, dip, rake, magnitude=6.0, source='unknown', - mtype='unknown'): +def fill_tensor_from_angles( + strike, dip, rake, magnitude=6.0, source="unknown", mtype="unknown" +): """Fill in moment tensor parameters from strike,dip, and rake. Args: @@ -37,13 +38,21 @@ def fill_tensor_from_angles(strike, dip, rake, magnitude=6.0, source='unknown', - rake (degrees) """ mt = plane_to_tensor(strike, dip, rake, mag=magnitude) - return fill_tensor_from_components(mt[0][0], mt[1][1], mt[2][2], - mt[1][0], mt[0][2], mt[1][2], - source=source, mtype=mtype) + return fill_tensor_from_components( + mt[0][0], + mt[1][1], + mt[2][2], + mt[1][0], + mt[0][2], + mt[1][2], + source=source, + mtype=mtype, + ) -def fill_tensor_from_components(mrr, mtt, mpp, mrt, mrp, mtp, source='unknown', - mtype='unknown'): +def fill_tensor_from_components( + mrr, mtt, mpp, mrt, mrp, mtp, source="unknown", mtype="unknown" +): """Fill in moment tensor parameters from moment tensor components. Args: @@ -56,7 +65,8 @@ def fill_tensor_from_components(mrr, mtt, mpp, mrt, mrp, mtp, source='unknown', source (str): Source (network, catalog) for input parameters. mtype (str): Focal mechanism or moment tensor type (Mww,Mwb,Mwc, etc.) Returns: - dict: Fully descriptive moment tensor dictionary, including fields: + tensor_params (dict): Fully descriptive moment tensor dictionary, including + fields: - mrr,mtt,mpp,mrt,mrp,mtp Moment tensor components. - T T-axis values: - azimuth (degrees) @@ -76,35 +86,27 @@ def fill_tensor_from_components(mrr, mtt, mpp, mrt, mrp, mtp, source='unknown', - dip (degrees) - rake (degrees) """ - tensor_params = {'mrr': mrr, - 'mtt': mtt, - 'mpp': mpp, - 'mrt': mrt, - 'mrp': mrp, - 'mtp': mtp} - tensor_params['source'] = source - tensor_params['type'] = mtype + tensor_params = { + "mrr": mrr, + "mtt": mtt, + "mpp": mpp, + "mrt": mrt, + "mrp": mrp, + "mtp": mtp, + } + tensor_params["source"] = source + tensor_params["type"] = mtype mt = MomentTensor(mrr, mtt, mpp, mrt, mrp, mtp, 1) tnp1 = mt2plane(mt) - np1 = {'strike': tnp1.strike, - 'dip': tnp1.dip, - 'rake': tnp1.rake} - tensor_params['NP1'] = np1.copy() - strike2, dip2, rake2 = aux_plane(np1['strike'], np1['dip'], np1['rake']) - np2 = {'strike': strike2, - 'dip': dip2, - 'rake': rake2} - tensor_params['NP2'] = np2.copy() + np1 = {"strike": tnp1.strike, "dip": tnp1.dip, "rake": tnp1.rake} + tensor_params["NP1"] = np1.copy() + strike2, dip2, rake2 = aux_plane(np1["strike"], np1["dip"], np1["rake"]) + np2 = {"strike": strike2, "dip": dip2, "rake": rake2} + tensor_params["NP2"] = np2.copy() T, N, P = mt2axes(mt) - tensor_params['T'] = {'azimuth': T.strike, - 'value': T.val, - 'plunge': T.dip} - tensor_params['N'] = {'azimuth': N.strike, - 'value': N.val, - 'plunge': N.dip} - tensor_params['P'] = {'azimuth': P.strike, - 'value': P.val, - 'plunge': P.dip} + tensor_params["T"] = {"azimuth": T.strike, "value": T.val, "plunge": T.dip} + tensor_params["N"] = {"azimuth": N.strike, "value": N.val, "plunge": N.dip} + tensor_params["P"] = {"azimuth": P.strike, "value": P.val, "plunge": P.dip} return tensor_params @@ -119,7 +121,7 @@ def plane_to_tensor(strike, dip, rake, mag=6.0): magnitude (float): Magnitude for moment tensor (not required if using moment tensor for angular comparisons.) Returns: - nparray: Tensor representation as 3x3 numpy matrix: + mt_matrix (nparray): Tensor representation as 3x3 numpy matrix: [[mrr, mrt, mrp] [mrt, mtt, mtp] [mrp, mtp, mpp]] @@ -132,22 +134,35 @@ def plane_to_tensor(strike, dip, rake, mag=6.0): # get tensor components mrr = mom * np.sin(2 * dip * d2r) * np.sin(rake * d2r) - mtt = -mom * ((np.sin(dip * d2r) * np.cos(rake * d2r) * np.sin(2 * strike * d2r)) + - (np.sin(2 * dip * d2r) * np.sin(rake * d2r) * - (np.sin(strike * d2r) * np.sin(strike * d2r)))) - mpp = mom * ((np.sin(dip * d2r) * np.cos(rake * d2r) * np.sin(2 * strike * d2r)) - - (np.sin(2 * dip * d2r) * np.sin(rake * d2r) * - (np.cos(strike * d2r) * np.cos(strike * d2r)))) - mrt = -mom * ((np.cos(dip * d2r) * np.cos(rake * d2r) * np.cos(strike * d2r)) + - (np.cos(2 * dip * d2r) * np.sin(rake * d2r) * np.sin(strike * d2r))) - mrp = mom * ((np.cos(dip * d2r) * np.cos(rake * d2r) * np.sin(strike * d2r)) - - (np.cos(2 * dip * d2r) * np.sin(rake * d2r) * np.cos(strike * d2r))) - mtp = -mom * ((np.sin(dip * d2r) * np.cos(rake * d2r) * np.cos(2 * strike * d2r)) + - (0.5 * np.sin(2 * dip * d2r) * np.sin(rake * d2r) * - np.sin(2 * strike * d2r))) + mtt = -mom * ( + (np.sin(dip * d2r) * np.cos(rake * d2r) * np.sin(2 * strike * d2r)) + + ( + np.sin(2 * dip * d2r) + * np.sin(rake * d2r) + * (np.sin(strike * d2r) * np.sin(strike * d2r)) + ) + ) + mpp = mom * ( + (np.sin(dip * d2r) * np.cos(rake * d2r) * np.sin(2 * strike * d2r)) + - ( + np.sin(2 * dip * d2r) + * np.sin(rake * d2r) + * (np.cos(strike * d2r) * np.cos(strike * d2r)) + ) + ) + mrt = -mom * ( + (np.cos(dip * d2r) * np.cos(rake * d2r) * np.cos(strike * d2r)) + + (np.cos(2 * dip * d2r) * np.sin(rake * d2r) * np.sin(strike * d2r)) + ) + mrp = mom * ( + (np.cos(dip * d2r) * np.cos(rake * d2r) * np.sin(strike * d2r)) + - (np.cos(2 * dip * d2r) * np.sin(rake * d2r) * np.cos(strike * d2r)) + ) + mtp = -mom * ( + (np.sin(dip * d2r) * np.cos(rake * d2r) * np.cos(2 * strike * d2r)) + + (0.5 * np.sin(2 * dip * d2r) * np.sin(rake * d2r) * np.sin(2 * strike * d2r)) + ) - mt_matrix = np.array([[mrr, mrt, mrp], - [mrt, mtt, mtp], - [mrp, mtp, mpp]]) + mt_matrix = np.array([[mrr, mrt, mrp], [mrt, mtt, mtp], [mrp, mtp, mpp]]) mt_matrix = mt_matrix * 1e-7 # convert from dyne-cm to N-m return mt_matrix