Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/model/geometry.py: 71%

1# https://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Representation of geometrical objects built from vertices, faces, outlines, 

8etc. 

9''' 

10 

11import numpy as num 

12 

13from pyrocko.guts import Object, List 

14from pyrocko.guts_array import Array 

15from pyrocko.table import Table, LocationRecipe 

16 

17from .event import Event 

18 

19from logging import getLogger 

20 

21 

22logger = getLogger('model.geometry') 

23 

24 

25def reduce_array_dims(array): 

26 ''' 

27 Support function to reduce output array ndims from table 

28 ''' 

29 

30 if array.shape[0] == 1 and array.ndim > 1: 

31 return num.squeeze(array, axis=0) 

32 else: 

33 return array 

34 

35 

36class Geometry(Object): 

37 ''' 

38 Spatial (and temporal) distributed properties of an event 

39 

40 The Geometry object allows to store properties of an event in a spatial 

41 (and temporal) distributed manner. For a set of planes ("faces"), 

42 characterized by their corner points ("vertices"), properties are stored. 

43 Also information on the outline of the source are stored. 

44 ''' 

45 

46 properties = Table.T( 

47 default=Table.D(), 

48 help='Properties that should be displayable on the surface of the' 

49 ' geometry. If 2d time dependency in column directions.', 

50 optional=True) 

51 

52 vertices = Table.T( 

53 default=Table.D(), 

54 help='Vertices of the mesh of the geometry. ' 

55 'Expected to be (lat,lon,north,east,depth) for each vertex.', 

56 optional=True) 

57 

58 faces = Table.T( 

59 default=Table.D(), 

60 help='Face integer indexes to the respective vertices. ' 

61 'Indexes belonging to one polygon have to be given row-wise.', 

62 optional=True) 

63 

64 outlines = List.T( 

65 Table.T(), 

66 default=[], 

67 help='List of vertices of the mesh of the outlines of the geometry. ' 

68 'Expected to be (lat,lon,north,east,depth) for each vertex' 

69 '(outline).', 

70 optional=True) 

71 

72 event = Event.T(default=Event.D()) 

73 

74 times = Array.T( 

75 shape=(None,), 

76 dtype='float64', 

77 help='1d vector of times [s] wrt. event time for which ' 

78 'properties have value. Must have constant delta t', 

79 optional=True) 

80 

81 def __init__(self, **kwargs): 

82 Object.__init__(self, **kwargs) 

83 self._ntimes = None 

84 

85 @property 

86 def deltat(self): 

87 ''' 

88 Sampling rate of properties (time difference) [s] 

89 ''' 

90 

91 if self.times.size > 2: 

92 return self.times[1] - self.times[0] 

93 else: 

94 return 1. 

95 

96 @property 

97 def ntimes(self): 

98 if self._ntimes is None: 

99 if self.times is not None: 

100 self._ntimes = len(self.times) 

101 else: 

102 return 0 

103 

104 return self._ntimes 

105 

106 def time2idx(self, time): 

107 if self.times.size > 2: 

108 idx = int(round((time - self.times.min() - ( 

109 self.deltat / 2.)) / self.deltat)) 

110 

111 if idx < 0: 

112 return 0 

113 elif idx > self.ntimes: 

114 return self.ntimes 

115 else: 

116 return idx 

117 else: 

118 return 0 

119 

120 def set_vertices(self, vertices): 

121 self.vertices.add_recipe(LocationRecipe()) 

122 self.vertices.add_col(( 

123 'c5', '', 

124 ('ref_lat', 'ref_lon', 'north_shift', 'east_shift', 'depth')), 

125 vertices) 

126 

127 def get_vertices(self, col='c5'): 

128 if self.vertices: 

129 return reduce_array_dims(self.vertices.get_col(col)) 

130 

131 def set_faces(self, faces): 

132 ncorners = faces.shape[1] 

133 sub_headers = tuple(['f{}'.format(i) for i in range(ncorners)]) 

134 self.faces.add_col(('faces', '', sub_headers), faces) 

135 

136 def get_faces(self, col='faces'): 

137 if self.faces: 

138 return reduce_array_dims(self.faces.get_col(col)) 

139 

140 def no_faces(self): 

141 return self.faces.get_nrows() 

142 

143 def setup(self, vertices, faces, outlines=None): 

144 self.set_vertices(vertices) 

145 self.set_faces(faces) 

146 

147 if outlines is not None: 

148 self.set_outlines(outlines) 

149 

150 def set_outlines(self, outlines): 

151 if outlines is not None: 

152 self.outlines = [] 

153 for outline in outlines: 

154 outl = Table() 

155 outl.add_recipe(LocationRecipe()) 

156 outl.add_col(( 

157 'c5', '', 

158 ('ref_lat', 'ref_lon', 'north_shift', 'east_shift', 

159 'depth')), 

160 outline) 

161 self.outlines.append(outl) 

162 

163 def get_outlines(self): 

164 if self.outlines: 

165 return self.outlines 

166 else: 

167 logger.warning('No outlines set!') 

168 

169 def add_property(self, name, values): 

170 if values.ndim == 1 or values.shape[1] == 1: 

171 self.properties.add_col(name, values.reshape(-1,)) 

172 elif (values.ndim == 2) and (self.times is not None): 

173 assert values.shape[1] == self.times.shape[0] 

174 self.properties.add_col(name, values) 

175 else: 

176 raise AttributeError( 

177 'Please give either 1D array or the associated times first.') 

178 

179 def get_property(self, name): 

180 return reduce_array_dims(self.properties.get_col(name)) 

181 

182 def has_property(self, name): 

183 return self.properties.has_col(name)