injection: adoptions to fault extraction

siria/dataset/faults.py

@@ -43,16 +43,30 @@ def extract_faults(
         radius=None,
         strike=None,
         dip=None,
-        rake=None):
+        rake=None,
+        restrict_to_radius=False):
 
     if not wesn:
         if lat is not None and lon is not None and radius is not None:
             wesn = pod.radius_to_region(lat, lon, radius)
 
     if wesn:
-        faults = [
-            f for f in faults
-            if pod.points_in_region(num.vstack((f.lat, f.lon)).T, wesn).any()]
+        faults_old = [f for f in faults]
+        faults = []
+
+        for i, f in enumerate(faults_old):
+            mask = pod.points_in_region(num.vstack((f.lat, f.lon)).T, wesn)
+
+            if not mask.any():
+                continue
+
+            if restrict_to_radius:
+                idx = num.arange(mask.shape[0])[mask]
+
+                f.lat = [f.lat[i] for i in idx]
+                f.lon = [f.lon[i] for i in idx]
+
+            faults.append(f)
 
     if strike:
         faults = [
@@ -76,12 +90,17 @@ def extract_faults(
 
 
 def load_faults_from_database(
-        database='GEM', exclude_depthcontours=True, **kwargs):
-    kwargs_extract = {
-        k: v for k, v in kwargs.items()
-        if k in extract_faults.__code__.co_varnames}
-
-    kwargs = {k: v for k, v in kwargs.items() if k not in kwargs_extract}
+        database='GEM',
+        exclude_depthcontours=True,
+        wesn=(),
+        lat=None,
+        lon=None,
+        radius=None,
+        strike=None,
+        dip=None,
+        rake=None,
+        restrict_to_radius=False,
+        **kwargs):
 
     database = database.upper()
     if database == 'DISS':
@@ -100,8 +119,16 @@ def load_faults_from_database(
             'No fault database {} found. Feel free to implement it'.format(
                 database))
 
-    if kwargs_extract:
-        faults = extract_faults(faults, **kwargs_extract)
+    faults = extract_faults(
+        faults,
+        wesn=wesn,
+        lat=lat,
+        lon=lon,
+        radius=radius,
+        strike=strike,
+        dip=dip,
+        rake=rake,
+        restrict_to_radius=restrict_to_radius)
 
     return faults
 

siria/si/invert.py

@@ -222,11 +222,13 @@ class SourceInjector(Object):
             list(samples[:, 1]),
             list(samples[:, 2]))
 
-    def _extract_from_faults(self, lat, lon, mt=None, magnitude=None):
+    def _extract_from_faults(
+            self, lat, lon, mt=None, magnitude=None, grond_config=None):
         '''
         Get faults from chosen fault database for given location.
         '''
         conf = self.config
+
         if conf.fault_radius_max is not None:
             radius = conf.fault_radius_max
         elif mt is not None:
@@ -241,14 +243,23 @@ class SourceInjector(Object):
                 for rake in (0., 90., 270.)])
             logger.info('Caluclated fault search radius of %f m using '
                         'magnitude', radius)
-        else:
-            raise ValueError('Can not find fault search radius!')
+
+        rg = grond_config.problem_config.ranges
+
+        radius_grond = num.sqrt(
+            num.max([rg['north_shift'].start, rg['north_shift'].stop])**2. +
+            num.max([rg['east_shift'].start, rg['east_shift'].stop])**2.)
+
+        radius += radius_grond
+        logger.info('Lateral model space extent of %f m added to fault '
+                    'search radius', radius_grond)
 
         fault_list = faults.load_faults_from_database(
             database=conf.fault_database,
             lat=lat,
             lon=lon,
-            radius=radius)
+            radius=radius,
+            restrict_to_radius=True)
 
         lats, lons, strikes, dips, rakes = [], [], [], [], []
 
@@ -323,7 +334,8 @@ class SourceInjector(Object):
             lat,
             lon,
             mt=event.moment_tensor,
-            magnitude=event.magnitude)
+            magnitude=event.magnitude,
+            grond_config=gronf)
 
         norths += list(fnorths)
         easts += list(feasts)

test/dataset/test_faults.py

@@ -0,0 +1,73 @@
+import unittest
+
+import numpy as num
+
+from pyrocko.model import Event
+from pyrocko import orthodrome as pod
+from siria.dataset import faults
+
+km = 1e3
+
+
+class DefaultValues(object):
+    event = Event.load('''--- !pf.Event
+lat: 42.74
+lon: 13.12
+time: '2016-10-30 06:40:23.900000095'
+depth: 10000.0
+name: 'Norcia_20161030_064018'
+magnitude: 6.633188496075327
+region: 'central Italy'
+catalog: 'INGV'
+moment_tensor: !pf.MomentTensor
+  mnn: 2.0e+18
+  mee: 8.0e+18
+  mdd: -1.0e+19
+  mne: 3.1e+18
+  mnd: -1.0e+18
+  med: 2.3e+18
+  strike1: 148.16380329988945
+  dip1: 40.625674815226844
+  rake1: -104.0702038392869
+  strike2: 346.4380332066921
+  dip2: 50.83320936298834
+  rake2: -78.21926297803395
+  moment: 9.99499874937461e+18
+  magnitude: 6.633188496075327
+''')
+
+    radius = 12879.
+
+
+defaults = DefaultValues()
+
+
+class TsumapsTest(unittest.TestCase):
+    def test_database_loaded(self):
+        fault_list = faults.load_faults_from_database(
+            database='EDSF',
+            lat=defaults.event.lat,
+            lon=defaults.event.lon,
+            radius=defaults.radius,
+            restrict_to_radius=True)
+
+        lats, lons = [], []
+
+        for f in fault_list:
+            if any([p is None for p in (
+                    f.lat, f.lon,
+                    f.strike_min, f.strike_max,
+                    f.dip_min, f.dip_max,
+                    f.rake_min, f.rake_max)]):
+
+                continue
+
+            lats += list(num.repeat(f.lat, 3))
+            lons += list(num.repeat(f.lon, 3))
+
+        norths, easts = pod.latlon_to_ne_numpy(
+            defaults.event.lat, defaults.event.lon,
+            num.array(lats), num.array(lons))
+
+        assert all(num.abs(norths) <= defaults.radius)
+        assert all(num.abs(easts) <= defaults.radius)