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=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)
46 xdists = dists[kept][:, kept]
47 xbadnesses = badnesses[kept]
48 xdeleted = _weed(xdists, xbadnesses, neighborhood, interaction_radius,
49 del_frac, max_del-ndeleted, max_depth, depth+1)
51 deleted[kept] = xdeleted
52 return deleted
55def weed(x, y, badnesses, neighborhood=1, nwanted=None, interaction_radius=3.):
56 assert x.size == y.size
57 n = x.size
58 NEW = num.newaxis
59 ax = x[NEW, :]
60 ay = y[NEW, :]
61 bx = x[:, NEW]
62 by = y[:, NEW]
64 if nwanted is None:
65 nwanted = n // 2
67 dx = num.abs(ax-bx)
68 dx = num.where(dx > 180., 360.-dx, dx)
70 dists = num.sqrt(dx**2 + (ay-by)**2)
71 deleted = _weed(
72 dists, badnesses, neighborhood, interaction_radius, del_frac=4,
73 max_del=n-nwanted, max_depth=500, depth=0)
75 kept = num.logical_not(deleted)
76 dists_kept = dists[kept][:, kept]
77 meandists_kept = neighborhood_density(dists_kept, neighborhood)
78 return deleted, meandists_kept
81def badnesses_c_mean(badnesses_nslc):
82 # convert stream badnesses to station badnesses by averaging
83 badnesses_nsl = {}
84 for nslc, b in badnesses_nslc.items():
85 nsl = nslc[:3]
86 badnesses_nsl.setdefault(nsl, []).append(b)
88 for nsl in list(badnesses_nsl.keys()):
89 this_station_badness = badnesses_nsl[nsl]
90 badnesses_nsl[nsl] = sum(this_station_badness) \
91 / len(this_station_badness)
93 return badnesses_nsl
96def weed_stations(stations, nwanted, neighborhood=3, default_badness=1.0,
97 badnesses=None,
98 badnesses_ns={}, badnesses_nsl={}, badnesses_nslc={}):
100 azimuths = num.zeros(len(stations), dtype=float)
101 dists = num.zeros(len(stations), dtype=float)
102 for ista, sta in enumerate(stations):
103 azimuths[ista] = sta.azimuth
104 dists[ista] = sta.dist_deg
106 badnesses_nslc_mean_c = badnesses_c_mean(badnesses_nslc)
107 badnesses2 = num.ones(len(stations), dtype=float)
108 for ista, sta in enumerate(stations):
109 nsl = sta.network, sta.station, sta.location
110 if badnesses is not None:
111 b = badnesses[ista]
112 else:
113 b = 1.0
114 b *= badnesses_ns.get(nsl, 1.0)
115 b *= badnesses_nsl.get(nsl, 1.0)
116 b *= badnesses_nslc_mean_c.get(nsl, 1.0)
117 badnesses2[ista] = b
119 deleted, meandists_kept = weed(
120 azimuths, dists, badnesses2,
121 nwanted=nwanted,
122 neighborhood=neighborhood)
124 stations_weeded = [station for (delete, station) in zip(deleted, stations)
125 if not delete]
127 return stations_weeded, meandists_kept, deleted