update according to pytomo3d

Author: wjlei1990

examples/adjoint_sources/parfile/multitaper.adjoint.50_100.config.yml

@@ -1,52 +1,53 @@
-# adjoint source type
-adj_src_type: "multitaper_misfit"
-
-# min and max period(unit: second)
-min_period: 50.0
-max_period: 100.0
-
-# adjoint config parameter
-lnpt: 15
-transfunc_waterlevel: 1.0E-10
-water_threshold: 0.02
-ipower_costaper: 10
-min_cycle_in_window: 3
-taper_percentage: 0.3
-mt_nw: 4.0
-num_taper: 5
-phase_step: 1.5
-dt_fac: 2.0
-err_fac: 2.5
-dt_max_scale: 3.5
-measure_type: 'dt'
-taper_type: 'hann'
-dt_sigma_min: 1.0
-dlna_sigma_min: 0.5
-use_cc_error: True
-use_mt_error: False
-
----
-# for postprocessing adjoint sources
-# interpolation. Starttime will be automatically set as cmt_time - 1.5 * hdur
-# to fit the SPECFEM behaviour
-interp_flag: True
-interp_delta: 0.1425
-interp_npts: 42000
-
-# for sum multiple insturments, like "II.AAK.00.BHZ" and "II.AAK.10.BHZ". if you turn
-# the weight_flag to be true, then you need also provide the weight_dict in the code
-sum_over_comp_flag: False
-weight_flag: True
-
-# filter the adjoint source
-filter_flag: True
-pre_filt: [0.0067, 0.01, 0.02, 0.025]
-
-# add missing components with zero trace(to prepare rotate)
-# If set to False, only rotate those with Both "R" and "T"
-# components. Set to True highly recommended unless you are
-# sreu what you are doing
-add_missing_comp_flag: False
-
-# rotate the adjoint source
-rotate_flag: False
+adjoint_config:
+  # adjoint source type
+  adj_src_type: "multitaper_misfit"
+  
+  # min and max period(unit: second)
+  min_period: 50.0
+  max_period: 100.0
+  
+  # adjoint config parameter
+  lnpt: 15
+  transfunc_waterlevel: 1.0E-10
+  water_threshold: 0.02
+  ipower_costaper: 10
+  min_cycle_in_window: 3
+  taper_percentage: 0.3
+  mt_nw: 4.0
+  num_taper: 5
+  phase_step: 1.5
+  dt_fac: 2.0
+  err_fac: 2.5
+  dt_max_scale: 3.5
+  measure_type: 'dt'
+  taper_type: 'hann'
+  dt_sigma_min: 1.0
+  dlna_sigma_min: 0.5
+  use_cc_error: True
+  use_mt_error: False
+  
+process_config:
+  # for postprocessing adjoint sources
+  # interpolation. Starttime will be automatically set as cmt_time - 1.5 * hdur
+  # to fit the SPECFEM behaviour
+  interp_flag: True
+  interp_delta: 0.1425
+  interp_npts: 42000
+  
+  # for sum multiple insturments, like "II.AAK.00.BHZ" and "II.AAK.10.BHZ". if you turn
+  # the weight_flag to be true, then you need also provide the weight_dict in the code
+  sum_over_comp_flag: False
+  weight_flag: True
+  
+  # filter the adjoint source
+  filter_flag: True
+  pre_filt: [0.0067, 0.01, 0.02, 0.025]
+  
+  # add missing components with zero trace(to prepare rotate)
+  # If set to False, only rotate those with Both "R" and "T"
+  # components. Set to True highly recommended unless you are
+  # sreu what you are doing
+  add_missing_comp_flag: False
+  
+  # rotate the adjoint source
+  rotate_flag: False

pypaw/bins/filter_windows.py

@@ -16,24 +16,11 @@
 from __future__ import (absolute_import, division, print_function)
 import os
 import argparse
-from pytomo3d.window import filter_windows
+from pytomo3d.window.filter_windows import filter_windows, count_windows
 from .utils import load_json, dump_json, load_yaml
 
 
-def main():
-
-    parser = argparse.ArgumentParser()
-    parser.add_argument('-f', action='store', dest='path_file', required=True,
-                        help="path file")
-    parser.add_argument('-p', action='store', dest='param_file', required=True,
-                        help="param file")
-    parser.add_argument('-v', action='store_true', dest='verbose',
-                        help="verbose flag")
-    args = parser.parse_args()
-
-    paths = load_json(args.path_file)
-    params = load_yaml(args.param_file)
-
+def run_window_filter(paths, params, verbose=False):
     window_file = paths["window_file"]
     station_file = paths["station_file"]
     output_file = paths["output_file"]
@@ -45,12 +32,15 @@ def main():
     print("output filtered window file: %s" % output_file)
 
     windows = load_json(window_file)
+    nchans_old, nwins_old = count_windows(windows)
     stations = load_json(station_file)
     measurements = load_json(measurement_file)
 
     # filter the window based on given sensor types
     windows_new, log = filter_windows(
-        windows, stations, measurements, params, verbose=args.verbose)
+        windows, stations, measurements, params, verbose=verbose)
+
+    nchans_new, nwins_new = count_windows(windows_new)
     # dump the new windows file to replace the original one
     dump_json(windows_new, output_file)
 
@@ -59,6 +49,27 @@ def main():
     print("Log file located at: %s" % logfile)
     dump_json(log, logfile)
 
+    print("=" * 10 + " Summary " + "=" * 10)
+    print("channels: %d --> %d" % (nchans_old, nchans_new))
+    print("windows: %d -- > %d" % (nwins_old, nwins_new))
+
+
+def main():
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-f', action='store', dest='path_file', required=True,
+                        help="path file")
+    parser.add_argument('-p', action='store', dest='param_file', required=True,
+                        help="param file")
+    parser.add_argument('-v', action='store_true', dest='verbose',
+                        help="verbose flag")
+    args = parser.parse_args()
+
+    paths = load_json(args.path_file)
+    params = load_yaml(args.param_file)
+
+    run_window_filter(paths, params, verbose=args.verbose)
+
 
 if __name__ == "__main__":
     main()

pypaw/sum_adjoint.py

@@ -103,10 +103,7 @@ def check_event_information_in_asdf_files(asdf_files):
         ds = ASDFDataSet(asdf_fn, mode='r')
         asdf_events[asdf_fn] = ds.events
 
-    diffs = check_events_consistent(asdf_events)
-    if len(diffs) != 0:
-        raise ValueError("Event information in %s not the same as others: %s"
-                         % (diffs, asdf_events.keys()))
+    check_events_consistent(asdf_events)
 
     event_base = asdf_events[asdf_events.keys()[0]]
     origin = event_base[0].preferred_origin()
@@ -164,6 +161,7 @@ def attach_station_to_db(self, station_info):
     def add_adjoint_dataset_on_channel_weight(self, ds, weights):
         """
         Add adjoint source based on channel window weight
+        The adjoint source's channel should be renamed to "MX"
         """
         misfits = {}
         adjsrc_group = ds.auxiliary_data.AdjointSources

pypaw/window_weights.py

@@ -20,9 +20,9 @@
 plt.switch_backend('agg')  # NOQA
 
 from pyasdf import ASDFDataSet
-from spaceweight import SpherePoint
-from spaceweight import SphereDistRel
 from pytomo3d.adjoint.sum_adjoint import check_events_consistent
+from pytomo3d.window.window_weights import determine_receiver_weighting, \
+    determine_category_weighting
 from pypaw.bins.utils import load_json, dump_json, load_yaml
 
 
@@ -147,125 +147,6 @@ def extract_source_location(input_info):
     return src_info
 
 
-def _receiver_validator(weights, rec_wcounts, cat_wcounts):
-
-    for comp in weights:
-        wsum = 0
-        for trace_id, trace_weight in weights[comp].iteritems():
-            nwin = rec_wcounts[trace_id]
-            wsum += trace_weight * nwin
-
-        if not np.isclose(wsum, cat_wcounts[comp]):
-            raise ValueError("receiver validator fails: %f, %f" %
-                             (wsum, cat_wcounts[comp]))
-
-
-def determine_receiver_weighting(src, stations, windows, max_ratio=0.35,
-                                 flag=True, plot=False,
-                                 figname_prefix=None):
-    """
-    Given one station and window information, determine the receiver
-    weighting
-    In one asdf file, there are still 3 components, for example,
-    ["BHR", "BHT", "BHZ"]. These three components should be treated
-    indepandently and weights will be calculated independantly.
-
-    :return: dict of weights which contains 3 components. Each components
-        contains weights values
-    """
-    center = SpherePoint(src["latitude"], src["longitude"],
-                         tag="source")
-
-    # extract window information
-    weights = {}
-    rec_wcounts = {}
-    cat_wcounts = defaultdict(lambda: 0)
-    for sta, sta_window in windows.iteritems():
-        for chan, chan_win in sta_window.iteritems():
-            comp = chan.split(".")[-1]
-            _nwin = len(chan_win)
-            if _nwin == 0:
-                continue
-            weights.setdefault(comp, {}).update({chan: 0.0})
-            rec_wcounts[chan] = _nwin
-            cat_wcounts[comp] += _nwin
-
-    # in each components, calculate weight
-    ref_dists = {}
-    cond_nums = {}
-    for comp, comp_info in weights.iteritems():
-        points = []
-        logger.info("Components:%s" % comp)
-        for chan in comp_info:
-            _comp = chan[-1]
-            if _comp == "Z":
-                point = SpherePoint(stations[chan]["latitude"],
-                                    stations[chan]["longitude"],
-                                    tag=chan, weight=1.0)
-            else:
-                # for R and T component. In station file, there
-                # are only E and N component. So we need transfer
-                echan = chan[:-1] + "E"
-                chan1 = chan[:-1] + "1"
-                zchan = chan[:-1] + "Z"
-                if echan in stations:
-                    point = SpherePoint(stations[echan]["latitude"],
-                                        stations[echan]["longitude"],
-                                        tag=chan, weight=1.0)
-                elif chan1 in stations:
-                    point = SpherePoint(stations[chan1]["latitude"],
-                                        stations[chan1]["longitude"],
-                                        tag=chan, weight=1.0)
-                elif zchan in stations:
-                    point = SpherePoint(stations[zchan]["latitude"],
-                                        stations[zchan]["longitude"],
-                                        tag=chan, weight=1.0)
-            points.append(point)
-
-        if flag:
-            # calculate weight; otherwise, leave it as default value(1)
-            weightobj = SphereDistRel(points, center=center)
-            scan_figname = figname_prefix + ".%s.smart_scan.png" % comp
-            ref_dists[comp], cond_nums[comp] = weightobj.smart_scan(
-                max_ratio=max_ratio, start=0.5, gap=0.5,
-                drop_ratio=0.95, plot=plot,
-                figname=scan_figname)
-            if plot:
-                figname = figname_prefix + ".%s.weight.png" % comp
-                weightobj.plot_global_map(figname=figname, lon0=180.0)
-        else:
-            ref_dists[comp] = None
-            cond_nums[comp] = None
-
-        wsum = 0
-        for point in points:
-            nwin = rec_wcounts[point.tag]
-            wsum += point.weight * nwin
-        norm_factor = cat_wcounts[comp] / wsum
-
-        for point in points:
-            weights[comp][point.tag] = point.weight * norm_factor
-
-    _receiver_validator(weights, rec_wcounts, cat_wcounts)
-
-    return {"rec_weights": weights, "rec_wcounts": rec_wcounts,
-            "cat_wcounts": cat_wcounts, "rec_ref_dists": ref_dists,
-            "rec_cond_nums": cond_nums}
-
-
-def _category_validator(weights, counts):
-    wsum = 0.0
-    nwins = 0
-    for p, pinfo in weights.iteritems():
-        for c in pinfo:
-            wsum += weights[p][c] * counts[p][c]
-            nwins += counts[p][c]
-
-    if not np.isclose(wsum, nwins):
-        raise ValueError("Category validator fails: %f, %f" %
-                         (wsum, nwins))
-
-
 def check_cat_consistency(cat_ratio, cat_wcounts):
     err = 0
     # check consistency
@@ -281,38 +162,6 @@ def check_cat_consistency(cat_ratio, cat_wcounts):
                          "consistent with window information")
 
 
-def determine_category_weighting(weight_param, cat_wcounts):
-    """
-    determine the category weighting based on window counts in each category
-    """
-    logger_block("Category Weighting")
-    weights = {}
-
-    cat_ratio = weight_param["ratio"]
-
-    print("cat_ratio: %s" % cat_ratio)
-    print("cat_wcounts: %s" % cat_wcounts)
-    sumv = 0
-    nwins = 0
-    for p, pinfo in cat_wcounts.iteritems():
-        for c in pinfo:
-            sumv += cat_wcounts[p][c] / cat_ratio[p][c]
-            nwins += cat_wcounts[p][c]
-
-    normc = nwins / sumv
-    logger.info("Total number of windows: %d" % nwins)
-
-    weights = {}
-    for p, pinfo in cat_wcounts.iteritems():
-        weights[p] = {}
-        for c in pinfo:
-            weights[p][c] = normc / cat_ratio[p][c]
-
-    logger.info("Category weights: %s" % weights)
-    _category_validator(weights, cat_wcounts)
-    return weights
-
-
 def plot_histogram(figname, array, nbins=50):
     # plot histogram of weights
     plt.hist(array, nbins)
@@ -529,9 +378,9 @@ def calculate_receiver_weights_asdf(self, period_info, weighting_param):
 
         _results = determine_receiver_weighting(
             self.src_info, station_info, window_info,
-            max_ratio=search_ratio,
-            flag=weight_flag,
-            plot=plot_flag, figname_prefix=figname_prefix)
+            search_ratio=search_ratio,
+            weight_flag=weight_flag,
+            plot_flag=plot_flag, figname_prefix=figname_prefix)
 
         return _results