1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
8from pyrocko.guts import Object, List
9from pyrocko.guts_array import Array
10from pyrocko.table import Table, LocationRecipe
12from .event import Event
14from logging import getLogger
17logger = getLogger('model.geometry')
20def reduce_array_dims(array):
21 '''
22 Support function to reduce output array ndims from table
23 '''
25 if array.shape[0] == 1 and array.ndim > 1:
26 return num.squeeze(array, axis=0)
27 else:
28 return array
31class Geometry(Object):
32 '''
33 Spatial (and temporal) distributed properties of an event
35 The Geometry object allows to store properties of an event in a spatial
36 (and temporal) distributed manner. For a set of planes ("faces"),
37 characterized by their corner points ("vertices"), properties are stored.
38 Also information on the outline of the source are stored.
39 '''
41 properties = Table.T(
42 default=Table.D(),
43 help='Properties that should be displayable on the surface of the'
44 ' geometry. If 2d time dependency in column directions.',
45 optional=True)
47 vertices = Table.T(
48 default=Table.D(),
49 help='Vertices of the mesh of the geometry. '
50 'Expected to be (lat,lon,north,east,depth) for each vertex.',
51 optional=True)
53 faces = Table.T(
54 default=Table.D(),
55 help='Face integer indexes to the respective vertices. '
56 'Indexes belonging to one polygon have to be given row-wise.',
57 optional=True)
59 outlines = List.T(
60 Table.T(),
61 default=[],
62 help='List of vertices of the mesh of the outlines of the geometry. '
63 'Expected to be (lat,lon,north,east,depth) for each vertex'
64 '(outline).',
65 optional=True)
67 event = Event.T(default=Event.D())
69 times = Array.T(
70 shape=(None,),
71 dtype='float64',
72 help='1d vector of times [s] wrt. event time for which '
73 'properties have value. Must have constant delta t',
74 optional=True)
76 def __init__(self, **kwargs):
77 Object.__init__(self, **kwargs)
78 self._ntimes = None
80 @property
81 def deltat(self):
82 '''
83 Sampling rate of properties (time difference) [s]
84 '''
86 if self.times.size > 2:
87 return self.times[1] - self.times[0]
88 else:
89 return 1.
91 @property
92 def ntimes(self):
93 if self._ntimes is None:
94 if self.times is not None:
95 self._ntimes = len(self.times)
96 else:
97 return 0
99 return self._ntimes
101 def time2idx(self, time):
102 if self.times.size > 2:
103 idx = int(round((time - self.times.min() - (
104 self.deltat / 2.)) / self.deltat))
106 if idx < 0:
107 return 0
108 elif idx > self.ntimes:
109 return self.ntimes
110 else:
111 return idx
112 else:
113 return 0
115 def set_vertices(self, vertices):
116 self.vertices.add_recipe(LocationRecipe())
117 self.vertices.add_col((
118 'c5', '',
119 ('ref_lat', 'ref_lon', 'north_shift', 'east_shift', 'depth')),
120 vertices)
122 def get_vertices(self, col='c5'):
123 if self.vertices:
124 return reduce_array_dims(self.vertices.get_col(col))
126 def set_faces(self, faces):
127 ncorners = faces.shape[1]
128 sub_headers = tuple(['f{}'.format(i) for i in range(ncorners)])
129 self.faces.add_col(('faces', '', sub_headers), faces)
131 def get_faces(self, col='faces'):
132 if self.faces:
133 return reduce_array_dims(self.faces.get_col(col))
135 def no_faces(self):
136 return self.faces.get_nrows()
138 def setup(self, vertices, faces, outlines=None):
139 self.set_vertices(vertices)
140 self.set_faces(faces)
142 if outlines is not None:
143 self.set_outlines(outlines)
145 def set_outlines(self, outlines):
146 if outlines is not None:
147 self.outlines = []
148 for outline in outlines:
149 outl = Table()
150 outl.add_recipe(LocationRecipe())
151 outl.add_col((
152 'c5', '',
153 ('ref_lat', 'ref_lon', 'north_shift', 'east_shift',
154 'depth')),
155 outline)
156 self.outlines.append(outl)
158 def get_outlines(self):
159 if self.outlines:
160 return self.outlines
161 else:
162 logger.warning('No outlines set!')
164 def add_property(self, name, values):
165 if values.ndim == 1 or values.shape[1] == 1:
166 self.properties.add_col(name, values.reshape(-1,))
167 elif (values.ndim == 2) and (self.times is not None):
168 assert values.shape[1] == self.times.shape[0]
169 self.properties.add_col(name, values)
170 else:
171 raise AttributeError(
172 'Please give either 1D array or the associated times first.')
174 def get_property(self, name):
175 return reduce_array_dims(self.properties.get_col(name))
177 def has_property(self, name):
178 return self.properties.has_col(name)