Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/scenario/base.py: 75%
167 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-04 09:52 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Base classes for the scenario generators.
8'''
10import math
11import hashlib
12import logging
14import numpy as num
16from pyrocko.guts import Object, Int, Bool, Float
17from pyrocko import orthodrome as od
18from pyrocko.dataset import gshhg, topo
19from .error import ScenarioError, LocationGenerationError
21logger = logging.getLogger('pyrocko.scenario.base')
23guts_prefix = 'pf.scenario'
25km = 1000.
26d2r = num.pi/180.
27N = 10000000
29coastlines = None
32def get_gsshg():
33 global coastlines
34 if coastlines is None:
35 logger.debug('Initialising GSHHG database.')
36 coastlines = gshhg.GSHHG.intermediate()
37 return coastlines
40def is_on_land(lat, lon, method='coastlines'):
41 if method == 'topo':
42 elevation = topo.elevation(lat, lon)
43 if elevation is None:
44 return False
45 else:
46 return topo.elevation(lat, lon) > 0.
48 elif method == 'coastlines':
49 is_land = get_gsshg().is_point_on_land(lat, lon)
50 logger.debug(
51 'Testing %.4f %.4f: %s' % (lat, lon, 'dry' if is_land else 'wet'))
53 return is_land
56def random_lat(rstate, lat_min=-90., lat_max=90.):
57 lat_min_ = 0.5*(math.sin(lat_min * math.pi/180.)+1.)
58 lat_max_ = 0.5*(math.sin(lat_max * math.pi/180.)+1.)
59 return math.asin(rstate.uniform(lat_min_, lat_max_)*2.-1.)*180./math.pi
62def random_latlon(rstate, avoid_water, ntries, label):
63 for itry in range(ntries):
64 logger.debug('%s: try %i' % (label, itry))
65 lat = random_lat(rstate)
66 lon = rstate.uniform(-180., 180.)
67 if not avoid_water or is_on_land(lat, lon):
68 return lat, lon
70 if avoid_water:
71 sadd = ' (avoiding water)'
73 raise LocationGenerationError('Could not generate location%s.' % sadd)
76def random_uniform(rstate, xmin, xmax, xdefault):
77 if None in (xmin, xmax):
78 return xdefault
79 else:
80 return rstate.uniform(xmin, xmax)
83class Generator(Object):
84 '''
85 Base class for random object generators in :py:mod:`pyrocko.scenario`.
86 '''
87 seed = Int.T(
88 optional=True,
89 help='Random seed for a reproducible scenario.')
91 def __init__(self, **kwargs):
92 Object.__init__(self, **kwargs)
93 self._seed = None
94 self._parent = None
95 self.update_hierarchy()
96 self._retry_offset = 0
98 def retry(self):
99 self.clear()
100 self._retry_offset += 1
101 for val in self.T.ivals(self):
102 if isinstance(val, Generator):
103 val.retry()
105 def clear(self):
106 self._seed = None
108 def hash(self):
109 return hashlib.sha1(
110 (self.dump() + '\n\n%i' % self._retry_offset).encode('utf8'))\
111 .hexdigest()
113 def get_seed_offset(self):
114 return int(self.hash(), base=16) % N
116 def update_hierarchy(self, parent=None):
117 self._parent = parent
118 for val in self.T.ivals(self):
119 if isinstance(val, Generator):
120 val.update_hierarchy(parent=self)
121 elif isinstance(val, list):
122 for el in val:
123 if isinstance(el, Generator):
124 el.update_hierarchy(parent=self)
126 def get_root(self):
127 if self._parent is None:
128 return self
129 else:
130 return self._parent.get_root()
132 def get_seed(self):
133 if self._seed is None:
134 if self.seed is None:
135 if self._parent is not None:
136 self._seed = self._parent.get_seed()
137 else:
138 self._seed = num.random.randint(N)
139 elif self.seed == 0:
140 self._seed = num.random.randint(N)
141 else:
142 self._seed = self.seed
144 return self._seed + self.get_seed_offset()
146 def get_rstate(self, i):
147 return num.random.RandomState(int(self.get_seed() + i))
149 def get_center_latlon(self):
150 return self._parent.get_center_latlon()
152 def get_radius(self):
153 return self._parent.get_radius()
155 def get_stations(self):
156 return []
159class LocationGenerator(Generator):
160 '''
161 Base class for generators providing random locations.
162 '''
164 avoid_water = Bool.T(
165 default=True,
166 help='Set whether wet areas should be avoided.')
167 center_lat = Float.T(
168 optional=True,
169 help='Center latitude for the random locations in [deg].')
170 center_lon = Float.T(
171 optional=True,
172 help='Center longitude for the random locations in [deg].')
173 radius = Float.T(
174 optional=True,
175 help='Radius for distribution of random locations [m].')
176 ntries = Int.T(
177 default=10,
178 help='Maximum number of tries to find a location satisifying all '
179 'given constraints')
180 north_shift_min = Float.T(
181 optional=True,
182 help='If given, lower bound of random northward cartesian offset [m].')
183 north_shift_max = Float.T(
184 optional=True,
185 help='If given, upper bound of random northward cartesian offset [m].')
186 east_shift_min = Float.T(
187 optional=True,
188 help='If given, lower bound of random eastward cartesian offset [m].')
189 east_shift_max = Float.T(
190 optional=True,
191 help='If given, upper bound of random eastward cartesian offset [m].')
192 depth_min = Float.T(
193 optional=True,
194 help='If given, minimum depth [m].')
195 depth_max = Float.T(
196 optional=True,
197 help='If given, maximum depth [m].')
199 def __init__(self, **kwargs):
200 Generator.__init__(self, **kwargs)
201 self._center_latlon = None
203 def clear(self):
204 Generator.clear(self)
205 self._center_latlon = None
207 def get_center_latlon(self):
208 if (self.center_lat is None) != (self.center_lon is None):
209 raise ScenarioError(
210 'Set both: lat and lon, or neither of them (in %s).'
211 % self.__class__.__name__)
213 if self._center_latlon is None:
215 if self.center_lat is not None and self.center_lon is not None:
216 self._center_latlon = self.center_lat, self.center_lon
218 else:
219 if self._parent:
220 self._center_latlon = self._parent.get_center_latlon()
221 else:
222 rstate = self.get_rstate(0)
223 self._center_latlon = random_latlon(
224 rstate, self.avoid_water, self.ntries,
225 self.__class__.__name__)
227 return self._center_latlon
229 def get_radius(self):
230 if self.radius is not None:
231 return self.radius
232 elif self._parent is not None:
233 return self._parent.get_radius()
234 else:
235 return None
237 def get_latlon(self, i):
238 rstate = self.get_rstate((i+1)*3+0)
239 for itry in range(self.ntries):
240 logger.debug('%s: try %i' % (self.__class__.__name__, itry))
241 radius = self.get_radius()
242 if radius is None:
243 lat = random_lat(rstate)
244 lon = rstate.uniform(-180., 180.)
245 elif radius == 0.0:
246 lat, lon = self.get_center_latlon()
247 else:
248 lat_center, lon_center = self.get_center_latlon()
249 while True:
250 north = rstate.uniform(-radius, radius)
251 east = rstate.uniform(-radius, radius)
252 if math.sqrt(north**2 + east**2) <= radius:
253 break
255 lat, lon = od.ne_to_latlon(lat_center, lon_center, north, east)
257 if not self.avoid_water or is_on_land(lat, lon):
258 logger.debug('location ok: %g, %g' % (lat, lon))
259 return lat, lon
261 if self.avoid_water:
262 sadd = ' (avoiding water)'
264 raise LocationGenerationError('Could not generate location%s.' % sadd)
266 def get_cartesian_offset(self, i):
267 rstate = self.get_rstate((i+1)*3+1)
268 north_shift = random_uniform(
269 rstate, self.north_shift_min, self.north_shift_max, 0.0)
270 east_shift = random_uniform(
271 rstate, self.east_shift_min, self.east_shift_max, 0.0)
273 return north_shift, east_shift
275 def get_depth(self, i):
276 rstate = self.get_rstate((i+1)*3+1)
277 return random_uniform(
278 rstate, self.depth_min, self.depth_max, 0.0)
280 def get_coordinates(self, i):
281 lat, lon = self.get_latlon(i)
282 north_shift, east_shift = self.get_cartesian_offset(i)
283 depth = self.get_depth(i)
284 return lat, lon, north_shift, east_shift, depth