1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5import numpy as num 

6 

7 

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 

13 

14 

15def _weed(dists, badnesses, neighborhood=1, interaction_radius=3., 

16 del_frac=4, max_del=100, max_depth=100, depth=0): 

17 

18 if depth > max_depth: 

19 assert False, 'max recursion depth reached' 

20 

21 meandists = neighborhood_density(dists, neighborhood) 

22 

23 order = meandists.argsort() 

24 candidates = order[:order.size//del_frac+1] 

25 badness_candidates = badnesses[candidates] 

26 order_badness = (-badness_candidates).argsort() 

27 

28 order[:order.size//del_frac+1] = \ 

29 order[:order.size//del_frac+1][order_badness] 

30 

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])): 

38 

39 deleted[ind] = True 

40 ndeleted += 1 

41 

42 if ndeleted == 0: 

43 return deleted 

44 

45 kept = num.logical_not(deleted).nonzero()[0] 

46 

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) 

51 

52 deleted[kept] = xdeleted 

53 return deleted 

54 

55 

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] 

64 

65 if nwanted is None: 

66 nwanted = n // 2 

67 

68 dx = num.abs(ax-bx) 

69 dx = num.where(dx > 180., 360.-dx, dx) 

70 

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) 

75 

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 

80 

81 

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) 

88 

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) 

93 

94 return badnesses_nsl 

95 

96 

97def weed_stations(stations, nwanted, neighborhood=3, default_badness=1.0, 

98 badnesses=None, 

99 badnesses_ns={}, badnesses_nsl={}, badnesses_nslc={}): 

100 

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 

106 

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 

119 

120 deleted, meandists_kept = weed( 

121 azimuths, dists, badnesses2, 

122 nwanted=nwanted, 

123 neighborhood=neighborhood) 

124 

125 stations_weeded = [station for (delete, station) in zip(deleted, stations) 

126 if not delete] 

127 

128 return stations_weeded, meandists_kept, deleted