1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5import numpy as num
8def neighborhood_density(dists, neighborhood=1):
9 sdists = dists.copy()
10 sdists.sort(axis=1)
11 meandists = num.mean(sdists[:, 1:1+neighborhood], axis=1)
12 return meandists
15def _weed(dists, badnesses, neighborhood=1, interaction_radius=3.,
16 del_frac=4, max_del=100, max_depth=100, depth=0):
18 if depth > max_depth:
19 assert False, 'max recursion depth reached'
21 meandists = neighborhood_density(dists, neighborhood)
23 order = meandists.argsort()
24 candidates = order[:order.size//del_frac+1]
25 badness_candidates = badnesses[candidates]
26 order_badness = (-badness_candidates).argsort()
28 order[:order.size//del_frac+1] = \
29 order[:order.size//del_frac+1][order_badness]
31 deleted = num.zeros(order.size, dtype=num.bool)
32 ndeleted = 0
33 for i, ind in enumerate(order):
34 if (i < order.size // del_frac // 2 + 1
35 and ndeleted < max_del
36 and num.all(
37 dists[ind, deleted] > interaction_radius*meandists[ind])):
39 deleted[ind] = True
40 ndeleted += 1
42 if ndeleted == 0:
43 return deleted
45 kept = num.logical_not(deleted).nonzero()[0]
47 xdists = dists[num.meshgrid(kept, kept)]
48 xbadnesses = badnesses[kept]
49 xdeleted = _weed(xdists, xbadnesses, neighborhood, interaction_radius,
50 del_frac, max_del-ndeleted, max_depth, depth+1)
52 deleted[kept] = xdeleted
53 return deleted
56def weed(x, y, badnesses, neighborhood=1, nwanted=None, interaction_radius=3.):
57 assert x.size == y.size
58 n = x.size
59 NEW = num.newaxis
60 ax = x[NEW, :]
61 ay = y[NEW, :]
62 bx = x[:, NEW]
63 by = y[:, NEW]
65 if nwanted is None:
66 nwanted = n // 2
68 dx = num.abs(ax-bx)
69 dx = num.where(dx > 180., 360.-dx, dx)
71 dists = num.sqrt(dx**2 + (ay-by)**2)
72 deleted = _weed(
73 dists, badnesses, neighborhood, interaction_radius, del_frac=4,
74 max_del=n-nwanted, max_depth=500, depth=0)
76 kept = num.logical_not(deleted).nonzero()[0]
77 dists_kept = dists[num.meshgrid(kept, kept)]
78 meandists_kept = neighborhood_density(dists_kept, neighborhood)
79 return deleted, meandists_kept
82def badnesses_c_mean(badnesses_nslc):
83 # convert stream badnesses to station badnesses by averaging
84 badnesses_nsl = {}
85 for nslc, b in badnesses_nslc.items():
86 nsl = nslc[:3]
87 badnesses_nsl.setdefault(nsl, []).append(b)
89 for nsl in list(badnesses_nsl.keys()):
90 this_station_badness = badnesses_nsl[nsl]
91 badnesses_nsl[nsl] = sum(this_station_badness) \
92 / len(this_station_badness)
94 return badnesses_nsl
97def weed_stations(stations, nwanted, neighborhood=3, default_badness=1.0,
98 badnesses=None,
99 badnesses_ns={}, badnesses_nsl={}, badnesses_nslc={}):
101 azimuths = num.zeros(len(stations), dtype=float)
102 dists = num.zeros(len(stations), dtype=float)
103 for ista, sta in enumerate(stations):
104 azimuths[ista] = sta.azimuth
105 dists[ista] = sta.dist_deg
107 badnesses_nslc_mean_c = badnesses_c_mean(badnesses_nslc)
108 badnesses2 = num.ones(len(stations), dtype=float)
109 for ista, sta in enumerate(stations):
110 nsl = sta.network, sta.station, sta.location
111 if badnesses is not None:
112 b = badnesses[ista]
113 else:
114 b = 1.0
115 b *= badnesses_ns.get(nsl, 1.0)
116 b *= badnesses_nsl.get(nsl, 1.0)
117 b *= badnesses_nslc_mean_c.get(nsl, 1.0)
118 badnesses2[ista] = b
120 deleted, meandists_kept = weed(
121 azimuths, dists, badnesses2,
122 nwanted=nwanted,
123 neighborhood=neighborhood)
125 stations_weeded = [station for (delete, station) in zip(deleted, stations)
126 if not delete]
128 return stations_weeded, meandists_kept, deleted