Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/weeding.py: 0%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7An algorithm to select station subsets with good coverage in azimuth and 

8distance. 

9''' 

10 

11import numpy as num 

12 

13 

14def neighborhood_density(dists, neighborhood=1): 

15 sdists = dists.copy() 

16 sdists.sort(axis=1) 

17 meandists = num.mean(sdists[:, 1:1+neighborhood], axis=1) 

18 return meandists 

19 

20 

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

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

23 

24 if depth > max_depth: 

25 assert False, 'max recursion depth reached' 

26 

27 meandists = neighborhood_density(dists, neighborhood) 

28 

29 order = meandists.argsort() 

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

31 badness_candidates = badnesses[candidates] 

32 order_badness = (-badness_candidates).argsort() 

33 

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

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

36 

37 deleted = num.zeros(order.size, dtype=bool) 

38 ndeleted = 0 

39 for i, ind in enumerate(order): 

40 if (i < order.size // del_frac // 2 + 1 

41 and ndeleted < max_del 

42 and num.all( 

43 dists[ind, deleted] > interaction_radius*meandists[ind])): 

44 

45 deleted[ind] = True 

46 ndeleted += 1 

47 

48 if ndeleted == 0: 

49 return deleted 

50 

51 kept = num.logical_not(deleted) 

52 xdists = dists[kept][:, kept] 

53 xbadnesses = badnesses[kept] 

54 xdeleted = _weed(xdists, xbadnesses, neighborhood, interaction_radius, 

55 del_frac, max_del-ndeleted, max_depth, depth+1) 

56 

57 deleted[kept] = xdeleted 

58 return deleted 

59 

60 

61def weed(x, y, badnesses, neighborhood=1, nwanted=None, interaction_radius=3.): 

62 assert x.size == y.size 

63 n = x.size 

64 NEW = num.newaxis 

65 ax = x[NEW, :] 

66 ay = y[NEW, :] 

67 bx = x[:, NEW] 

68 by = y[:, NEW] 

69 

70 if nwanted is None: 

71 nwanted = n // 2 

72 

73 dx = num.abs(ax-bx) 

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

75 

76 dists = num.sqrt(dx**2 + (ay-by)**2) 

77 deleted = _weed( 

78 dists, badnesses, neighborhood, interaction_radius, del_frac=4, 

79 max_del=n-nwanted, max_depth=500, depth=0) 

80 

81 kept = num.logical_not(deleted) 

82 dists_kept = dists[kept][:, kept] 

83 meandists_kept = neighborhood_density(dists_kept, neighborhood) 

84 return deleted, meandists_kept 

85 

86 

87def badnesses_c_mean(badnesses_nslc): 

88 # convert stream badnesses to station badnesses by averaging 

89 badnesses_nsl = {} 

90 for nslc, b in badnesses_nslc.items(): 

91 nsl = nslc[:3] 

92 badnesses_nsl.setdefault(nsl, []).append(b) 

93 

94 for nsl in list(badnesses_nsl.keys()): 

95 this_station_badness = badnesses_nsl[nsl] 

96 badnesses_nsl[nsl] = sum(this_station_badness) \ 

97 / len(this_station_badness) 

98 

99 return badnesses_nsl 

100 

101 

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

103 badnesses=None, 

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

105 

106 azimuths = num.zeros(len(stations), dtype=float) 

107 dists = num.zeros(len(stations), dtype=float) 

108 for ista, sta in enumerate(stations): 

109 azimuths[ista] = sta.azimuth 

110 dists[ista] = sta.dist_deg 

111 

112 badnesses_nslc_mean_c = badnesses_c_mean(badnesses_nslc) 

113 badnesses2 = num.ones(len(stations), dtype=float) 

114 for ista, sta in enumerate(stations): 

115 nsl = sta.network, sta.station, sta.location 

116 if badnesses is not None: 

117 b = badnesses[ista] 

118 else: 

119 b = 1.0 

120 b *= badnesses_ns.get(nsl, 1.0) 

121 b *= badnesses_nsl.get(nsl, 1.0) 

122 b *= badnesses_nslc_mean_c.get(nsl, 1.0) 

123 badnesses2[ista] = b 

124 

125 deleted, meandists_kept = weed( 

126 azimuths, dists, badnesses2, 

127 nwanted=nwanted, 

128 neighborhood=neighborhood) 

129 

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

131 if not delete] 

132 

133 return stations_weeded, meandists_kept, deleted