wip

.vscode/settings.json

@@ -0,0 +1,11 @@
+{
+    "python.testing.pytestArgs": [
+        "tests"
+    ],
+    "python.testing.unittestEnabled": false,
+    "python.testing.nosetestsEnabled": false,
+    "python.testing.pytestEnabled": true,
+    "python.pythonPath": "/usr/bin/python3",
+    "python.linting.flake8Enabled": true,
+    "python.linting.enabled": true
+}

.vscode/tasks.json

@@ -0,0 +1,13 @@
+{
+    // See https://go.microsoft.com/fwlink/?LinkId=733558
+    // for the documentation about the tasks.json format
+    "version": "2.0.0",
+    "tasks": [
+        {
+            "label": "install pip",
+            "type": "process",
+            "command": "${config:python.pythonPath}",
+            "args": ["-m", "pip", "install", "."]
+        }
+    ]
+}

pyproject.toml

@@ -0,0 +1,6 @@
+[build-system]
+requires = [
+    "setuptools>=42",
+    "wheel"
+]
+build-backend = "setuptools.build_meta"

setup.py

@@ -4,7 +4,7 @@ with open("README.md", "r", encoding="utf-8") as fh:
     long_description = fh.read()
 
 setuptools.setup(
-    name="pyrocko-das",
+    name="pyrocko_das",
     version="0.0.1",
     author="Marius Paul Isken",
     author_email="mi@gfz-potsdam.de",
@@ -20,7 +20,7 @@ setuptools.setup(
         "License :: OSI Approved :: MIT License",
         "Operating System :: OS Independent",
     ],
-    package_dir={"": "src"},
-    packages=setuptools.find_packages(where="src"),
+    package_dir={"pyrocko_das": "src"},
+    packages=["pyrocko_das"],
     python_requires=">=3.6",
 )

src/__init__.py

@@ -1,27 +0,0 @@
-import pyrocko
-import numpy as num
-
-from .utils import traces_to_numpy_and_meta
-
-
-def remove_idas_instrument_vibration(traces, copy=True, return_ref=False):
-    data, meta = traces_to_numpy_and_meta(traces)
-
-    ntraces_ref = int((meta.start_distance * 0.9) / meta.spatial_resolution)
-
-    ref = num.mean(data[:ntraces_ref], axis=1)
-
-    out_traces = []
-    for tr in traces:
-        if copy:
-            tr = tr.copy()
-        tr.ydata -= ref
-        out_traces.append(tr)
-
-    if return_ref:
-        trace_ref = tr.copy(data=False)
-        trace_ref.ydata = ref
-        trace_ref.station = "ref"
-        out_traces.append(trace_ref)
-
-    return out_traces

src/gf.py

@@ -0,0 +1,264 @@
+import numpy as num
+from scipy import interpolate
+from collections.abc import Iterable
+
+from pyrocko.guts import Object, Float, List, Tuple, String, Timestamp
+from pyrocko import gf
+from pyrocko.model import Location
+import pyrocko.orthodrome as od
+
+r2d = 180./num.pi
+km = 1e3
+
+META = {
+    'measure_length': None,
+    'start_distance': None,
+    'stop_distance': None,
+    'gauge_length': None,
+    'spatial_resolution': None,
+
+    'geo_lat': None,
+    'geo_lon': None,
+    'geo_elevation': None,
+
+    'channel': None,
+    'unit': None
+}
+
+
+class QuantityType(gf.meta.QuantityType):
+    choices = gf.meta.QuantityType.choices + ['strain', 'strain_rate']
+
+
+class Fiber(Object):
+    lat = Float.T(
+        default=0.)
+    lon = Float.T(
+        default=0.)
+
+    codes = Tuple.T(
+        4, String.T(), default=('', '', '', 'HSF'),
+        help='network, station, location and channel codes to be set on '
+             'the response trace. If station code is empty it will be filled'
+             ' by the channel number')
+    quantity = QuantityType.T(
+        default='strain',
+        help='Measurement quantity type. If not given, it is guessed from the '
+             'channel code. For some common cases, derivatives of the stored '
+             'quantities are supported by using finite difference '
+             'approximations (e.g. displacement to velocity or acceleration). '
+             '4th order central FD schemes are used.')
+
+    elevation = Float.T(
+        default=0.0,
+        help='station surface elevation in [m]')
+
+    store_id = gf.meta.StringID.T(
+        optional=True,
+        help='ID of Green\'s function store to use for the computation. '
+             'If not given, the processor may use a system default.')
+
+    sample_rate = Float.T(
+        optional=True,
+        help='sample rate to produce. '
+             'If not given the GF store\'s default sample rate is used. '
+             'GF store specific restrictions may apply.')
+
+    tmin = Timestamp.T(
+        optional=True,
+        help='time of first sample to request in [s]. '
+             'If not given, it is determined from the Green\'s functions.')
+    tmax = Timestamp.T(
+        optional=True,
+        help='time of last sample to request in [s]. '
+             'If not given, it is determined from the Green\'s functions.')
+    interpolation = gf.meta.InterpolationMethod.T(
+        default='multilinear',
+        help='Interpolation method between Green\'s functions. Supported are'
+                ' ``nearest_neighbor`` and ``multilinear``')
+
+    coordinates = List.T(
+        help='coordinates of the cable as ``pyrocko.model.Location`` or a tuple'
+             ' of (north_shift, east_shift, [elevation]).')
+
+    channel_spacing = Float.T(
+        default=4.,
+        help='Channel spacing [m].')
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+        self._locations = []
+        for coord in self.coordinates:
+            if isinstance(coord, Location):
+                loc = coord
+            elif isinstance(coord, Iterable) and len(coord) in (2, 3):
+                loc = Location(
+                    lat=self.lat,
+                    lon=self.lon,
+                    north_shift=coord[0],
+                    east_shift=coord[1],
+                    elevation=self.elevation if len(coord) == 2 else coord[2])
+            else:
+                raise AttributeError(
+                    'coordinates have to be a list of pyrocko.model.Location'
+                    ' or tuples of (north_shift, east_shifts, [elevation]),'
+                    ' not %s' % type(coord))
+
+            self._locations.append(loc)
+        assert len(self._locations) > 1, 'Fiber needs more than 1 coordinates'
+
+
+    @classmethod
+    def from_stations(cls, stations, **kwargs):
+        return cls(
+            coordinates=stations,
+            **kwargs)
+
+    @property
+    def distances(self):
+        locations = self._locations
+        nlocs = len(locations)
+
+        dists = [0.]
+        for il in range(1, nlocs):
+            dists.append(
+                locations[il-1].distance_3d_to(locations[il]))
+
+        return num.array(dists)
+
+    @property
+    def length(self):
+        return num.sum(self.distances)
+
+    @property
+    def nchannels(self):
+        dists = self.length
+        return int(dists // self.channel_spacing)
+
+    def interpolate_channels(self, start_channel=0):
+        dists = num.cumsum(self.distances)
+
+        interp_lat = interpolate.interp1d(
+            dists,
+            tuple(loc.effective_lat for loc in self._locations))
+        interp_lon = interpolate.interp1d(
+            dists,
+            tuple(loc.effective_lon for loc in self._locations))
+        interp_elevation = interpolate.interp1d(
+            dists,
+            tuple(loc.elevation for loc in self._locations))
+
+        interp_distances = num.arange(
+            0., dists.max(),
+            step=self.channel_spacing)
+
+        channels = num.arange(interp_distances.size) + start_channel
+        nchannels = channels.size
+
+        lats = interp_lat(interp_distances)
+        lons = interp_lon(interp_distances)
+        elevations = interp_elevation(interp_distances)
+
+        azis = num.empty(nchannels)
+        azis[:-1] = od.azimuth_numpy(
+            lats[:-1], lons[:-1],
+            lats[1:], lons[1:],
+        )
+        azis[-1] = azis[-2]
+
+        ds = num.full(nchannels, self.channel_spacing)
+        dips = -num.arctan2(
+            num.gradient(elevations), ds) * r2d
+
+        return lats, lons, elevations, azis, dips, channels
+
+    def get_targets(self):
+        lats, lons, elevations, azis, dips, channels = \
+            self.interpolate_channels()
+        nchannels = self.nchannels
+
+        quantity = self.quantity
+        if self.quantity in ('strain', 'strain_rate'):
+            quantity = 'displacement'
+
+        targets = []
+        for icha in range(nchannels):
+            codes = list(self.codes)
+            if not codes[1]:
+                codes[1] = '%05d' % icha
+
+            t = gf.Target(
+                lat=lats[icha],
+                lon=lons[icha],
+                elevation=elevations[icha],
+                azimuth=azis[icha],
+                dip=dips[icha],
+                north_shift=0.,
+                east_shift=0.,
+                sample_rate=self.sample_rate,
+
+                codes=codes,
+                store_id=self.store_id,
+                quantity=quantity,
+                interpolation=self.interpolation
+            )
+            targets.append(t)
+        return targets
+
+    def get_locations(self):
+        return self._locations
+
+
+class LocalEngine(gf.LocalEngine):
+
+    def process_fiber(self, source, fiber, **kwargs):
+        assert isinstance(source, Iterable), \
+            'Currently only one source is supported!'
+
+        targets = fiber.get_targets()
+        resp = self.process(source, targets, **kwargs)
+
+        traces = resp.pyrocko_traces()
+        ntraces = len(traces)
+
+        all_times = [(tr.tmin, tr.tmax) for tr in traces]
+        all_tmin = num.min(all_times)
+        all_tmax = num.max(all_times)
+
+        meta = {
+            'measure_length': fiber.length,
+            'start_distance': 0.,
+            'stop_distance': fiber.length,
+            'gauge_length': fiber.channel_spacing,
+            'spatial_resolution': fiber.channel_spacing
+        }
+
+        for icha, (tr, target) in enumerate(zip(traces, targets)):
+            tr.extend(tmin=all_tmin, tmax=all_tmax, fillmethod='repeat')
+            tr.meta = meta.copy()
+            tr.meta['channel'] = icha
+            tr.meta['geo_lat'] = target.lat
+            tr.meta['geo_lon'] = target.lon
+            tr.meta['geo_elevation'] = target.elevation
+
+        nsamples = set([tr.ydata.size for tr in traces])
+        assert len(nsamples) == 1
+        nsamples = nsamples.pop()
+
+        if fiber.quantity in ('strain', 'strain_rate'):
+            assert ntraces > 3, \
+                'Too few channels in fiber for finite-difference derivation'
+
+            trs_strain_grad = num.gradient(
+                num.array([tr.ydata for tr in traces]), fiber.channel_spacing,
+                axis=0)
+
+            for itr, tr in enumerate(traces):
+                tr.set_ydata(trs_strain_grad[itr])
+
+        if fiber.quantity == 'strain_rate':
+            for tr in traces:
+                tr.differentiate()
+
+        return traces

src/restitution.py

@@ -2,21 +2,26 @@ import numpy as num
 from .utils import traces_to_numpy_and_meta
 
 
-def de_to_e(traces, copy=False):
+def rad_to_de(traces, copy=False, m_per_rad=11.6e-9):
     out_traces = []
     for tr in traces:
+        assert 'gauge_length' in tr.meta, 'gauge_length not in metadata'
         if copy:
             tr = tr.copy()
-        tr.set_ydata(num.cumsum(tr.ydata) * tr.deltat)
+        tr.ydata *= m_per_rad
+        tr.ydata /= tr.meta['gauge_length'] * tr.deltat
+        tr.meta['unit'] = 'strain rate (m/s/m)'
         out_traces.append(tr)
 
-    return out_traces
-
 
-def de_to_velocity_static_slowness(traces, slowness, copy=False):
-    out_traces = de_to_e(traces, copy)
-    for tr in out_traces:
-        tr.set_ydata(tr.ydata / slowness)
+def de_to_e(traces, copy=False):
+    out_traces = []
+    for tr in traces:
+        if copy:
+            tr = tr.copy()
+        tr.set_ydata(num.cumsum(tr.ydata) / tr.deltat)
+        tr.meta['unit'] = 'strain (m/m)'
+        out_traces.append(tr)
 
     return out_traces
 
@@ -27,6 +32,17 @@ def de_to_acceleration_static_slowness(traces, slowness, copy=False):
         if copy:
             tr = tr.copy()
         tr.set_ydata(tr.ydata / slowness)
+        tr.meta['unit'] = 'acceleration (m/s^2)'
+        out_traces.append(tr)
+
+    return out_traces
+
+
+def de_to_velocity_static_slowness(traces, slowness, copy=False):
+    out_traces = de_to_e(traces, copy)
+    for tr in out_traces:
+        tr.set_ydata(tr.ydata / slowness)
+        tr.meta['unit'] = 'velocity (m/s)'
 
     return out_traces
 
@@ -38,6 +54,7 @@ def de_to_relative_displacement(traces, copy=False):
     data = num.cumsum(data, axis=0) * meta.spatial_resolution
     for itr, tr in enumerate(traces):
         tr.set_ydata(data[itr])
+        tr.meta['unit'] = 'displacement (m)'
 
     return traces
 
@@ -49,5 +66,6 @@ def de_to_relative_velocity(traces, copy=False):
     data = num.diff(data, n=1, axis=0) / meta.spatial_resolution
     for itr, tr in enumerate(traces):
         tr.set_ydata(data[itr])
+        tr.meta['unit'] = 'velocity (m/s)'
 
     return traces

tests/test_gf.py

@@ -0,0 +1,110 @@
+from pyrocko_das import gf
+from pyrocko.model import Location
+from pyrocko import trace
+import numpy as num
+
+km = 1e3
+
+
+def test_fiber():
+
+    fiber = gf.Fiber(
+        lat=0.5,
+        lon=0.5,
+        store_id='test',
+
+        coordinates=(
+            (0, 0),
+            (100, 0)
+        ),
+        channel_spacing=1.
+    )
+    assert fiber.distances.size == 2
+    assert fiber.distances[0] == 0.
+    assert fiber.distances[1] == 100.
+
+    assert fiber.length == 100.
+    assert fiber.nchannels == 100
+
+    northings, eastings, elevations, azis, dips, channels = \
+        fiber.interpolate_channels()
+    assert channels.size == 100.
+
+    num.testing.assert_equal(azis, 0.)
+    num.testing.assert_equal(dips, 0.)
+
+    targets = fiber.get_targets()
+
+    assert len(targets) == 100
+    assert len(targets) == fiber.nchannels
+
+    fiber = gf.Fiber(
+        lat=0.,
+        lon=0.,
+        store_id='test',
+
+        coordinates=(
+            (0, 0, 0.),
+            (-20, 0, 20.),
+            (-40, 0, 40.),
+        ),
+        channel_spacing=1.,
+        interpolation='multilinear'
+    )
+
+    northings, eastings, elevations, azis, dips, channels = \
+        fiber.interpolate_channels()
+    num.testing.assert_equal(azis, 180.)
+    num.testing.assert_equal(dips, -45.)
+
+    fiber = gf.Fiber(
+        lat=0.,
+        lon=0.,
+        store_id='test',
+
+        coordinates=(
+            Location(
+                north_shift=0.,
+                east_shift=0,
+                elevation=0),
+            Location(
+                north_shift=-100.,
+                east_shift=0,
+                elevation=100),
+        ),
+        channel_spacing=1.,
+        interpolation='multilinear'
+    )
+
+    northings, eastings, elevations, azis, dips, channels = \
+        fiber.interpolate_channels()
+    num.testing.assert_equal(azis, 180.)
+    num.testing.assert_equal(dips, -45.)
+
+
+def test_process_fiber():
+    import pyrocko.gf as pgf
+
+    fiber = gf.Fiber(
+        lat=0.0,
+        lon=0.0,
+        quantity='acceleration',
+        store_id='das_test',
+
+        coordinates=(
+            (1100, 0),
+            (2000, 1000)
+        ),
+        channel_spacing=4.
+    )
+    engine = gf.LocalEngine(use_config=True)
+
+    source = pgf.DCSource(
+        lat=0., lon=0.,
+        depth=2*km,
+        strike=45.,
+        dip=30.
+    )
+
+    traces = engine.process_fiber(source, fiber)
+    trace.snuffle(traces)

tests/test_remove_idas_vibration.py

@@ -1,9 +0,0 @@
-from pyrocko import io
-from pyrocko_das import remove_idas_instrument_vibration
-
-
-traces = io.read("")
-
-
-def test_remove_idas_vibration():
-    remove_idas_instrument_vibration(traces)

tests/test_restitution.py

@@ -0,0 +1,14 @@
+import numpy as num
+
+from pyrocko import io
+from pyrocko import trace
+
+
+def generate_data(ncha=100, deltat=0.001, cha_spacing=1.):
+    meta = io.tdms_idas.META_KEYS.copy()
+
+    cha_dist = num.arange(ncha) * cha_spacing
+    traces = []
+
+    for icha in range(ncha):
+        dist = cha_dist[icha]