1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import hashlib
8import logging
10import numpy as num
12from pyrocko.guts import Object, Int, Bool, Float
13from pyrocko import orthodrome as od
14from pyrocko.dataset import gshhg, topo
15from .error import ScenarioError, LocationGenerationError
17logger = logging.getLogger('pyrocko.scenario.base')
19guts_prefix = 'pf.scenario'
21km = 1000.
22d2r = num.pi/180.
23N = 10000000
25coastlines = None
28def get_gsshg():
29 global coastlines
30 if coastlines is None:
31 logger.debug('Initialising GSHHG database.')
32 coastlines = gshhg.GSHHG.intermediate()
33 return coastlines
36def is_on_land(lat, lon, method='coastlines'):
37 if method == 'topo':
38 elevation = topo.elevation(lat, lon)
39 if elevation is None:
40 return False
41 else:
42 return topo.elevation(lat, lon) > 0.
44 elif method == 'coastlines':
45 is_land = get_gsshg().is_point_on_land(lat, lon)
46 logger.debug(
47 'Testing %.4f %.4f: %s' % (lat, lon, 'dry' if is_land else 'wet'))
49 return is_land
52def random_lat(rstate, lat_min=-90., lat_max=90.):
53 lat_min_ = 0.5*(math.sin(lat_min * math.pi/180.)+1.)
54 lat_max_ = 0.5*(math.sin(lat_max * math.pi/180.)+1.)
55 return math.asin(rstate.uniform(lat_min_, lat_max_)*2.-1.)*180./math.pi
58def random_latlon(rstate, avoid_water, ntries, label):
59 for itry in range(ntries):
60 logger.debug('%s: try %i' % (label, itry))
61 lat = random_lat(rstate)
62 lon = rstate.uniform(-180., 180.)
63 if not avoid_water or is_on_land(lat, lon):
64 return lat, lon
66 if avoid_water:
67 sadd = ' (avoiding water)'
69 raise LocationGenerationError('Could not generate location%s.' % sadd)
72def random_uniform(rstate, xmin, xmax, xdefault):
73 if None in (xmin, xmax):
74 return xdefault
75 else:
76 return rstate.uniform(xmin, xmax)
79class Generator(Object):
80 seed = Int.T(
81 optional=True,
82 help='Random seed for a reproducible scenario.')
84 def __init__(self, **kwargs):
85 Object.__init__(self, **kwargs)
86 self._seed = None
87 self._parent = None
88 self.update_hierarchy()
89 self._retry_offset = 0
91 def retry(self):
92 self.clear()
93 self._retry_offset += 1
94 for val in self.T.ivals(self):
95 if isinstance(val, Generator):
96 val.retry()
98 def clear(self):
99 self._seed = None
101 def hash(self):
102 return hashlib.sha1(
103 (self.dump() + '\n\n%i' % self._retry_offset).encode('utf8'))\
104 .hexdigest()
106 def get_seed_offset(self):
107 return int(self.hash(), base=16) % N
109 def update_hierarchy(self, parent=None):
110 self._parent = parent
111 for val in self.T.ivals(self):
112 if isinstance(val, Generator):
113 val.update_hierarchy(parent=self)
114 elif isinstance(val, list):
115 for el in val:
116 if isinstance(el, Generator):
117 el.update_hierarchy(parent=self)
119 def get_root(self):
120 if self._parent is None:
121 return self
122 else:
123 return self._parent.get_root()
125 def get_seed(self):
126 if self._seed is None:
127 if self.seed is None:
128 if self._parent is not None:
129 self._seed = self._parent.get_seed()
130 else:
131 self._seed = num.random.randint(N)
132 elif self.seed == 0:
133 self._seed = num.random.randint(N)
134 else:
135 self._seed = self.seed
137 return self._seed + self.get_seed_offset()
139 def get_rstate(self, i):
140 return num.random.RandomState(int(self.get_seed() + i))
142 def get_center_latlon(self):
143 return self._parent.get_center_latlon()
145 def get_radius(self):
146 return self._parent.get_radius()
148 def get_stations(self):
149 return []
152class LocationGenerator(Generator):
154 avoid_water = Bool.T(
155 default=True,
156 help='Set whether wet areas should be avoided.')
157 center_lat = Float.T(
158 optional=True,
159 help='Center latitude for the random locations in [deg].')
160 center_lon = Float.T(
161 optional=True,
162 help='Center longitude for the random locations in [deg].')
163 radius = Float.T(
164 optional=True,
165 help='Radius for distribution of random locations [m].')
166 ntries = Int.T(
167 default=10,
168 help='Maximum number of tries to find a location satisifying all '
169 'given constraints')
170 north_shift_min = Float.T(
171 optional=True,
172 help='If given, lower bound of random northward cartesian offset [m].')
173 north_shift_max = Float.T(
174 optional=True,
175 help='If given, upper bound of random northward cartesian offset [m].')
176 east_shift_min = Float.T(
177 optional=True,
178 help='If given, lower bound of random eastward cartesian offset [m].')
179 east_shift_max = Float.T(
180 optional=True,
181 help='If given, upper bound of random eastward cartesian offset [m].')
182 depth_min = Float.T(
183 optional=True,
184 help='If given, minimum depth [m].')
185 depth_max = Float.T(
186 optional=True,
187 help='If given, maximum depth [m].')
189 def __init__(self, **kwargs):
190 Generator.__init__(self, **kwargs)
191 self._center_latlon = None
193 def clear(self):
194 Generator.clear(self)
195 self._center_latlon = None
197 def get_center_latlon(self):
198 if (self.center_lat is None) != (self.center_lon is None):
199 raise ScenarioError(
200 'Set both: lat and lon, or neither of them (in %s).'
201 % self.__class__.__name__)
203 if self._center_latlon is None:
205 if self.center_lat is not None and self.center_lon is not None:
206 self._center_latlon = self.center_lat, self.center_lon
208 else:
209 if self._parent:
210 self._center_latlon = self._parent.get_center_latlon()
211 else:
212 rstate = self.get_rstate(0)
213 self._center_latlon = random_latlon(
214 rstate, self.avoid_water, self.ntries,
215 self.__class__.__name__)
217 return self._center_latlon
219 def get_radius(self):
220 if self.radius is not None:
221 return self.radius
222 elif self._parent is not None:
223 return self._parent.get_radius()
224 else:
225 return None
227 def get_latlon(self, i):
228 rstate = self.get_rstate((i+1)*3+0)
229 for itry in range(self.ntries):
230 logger.debug('%s: try %i' % (self.__class__.__name__, itry))
231 radius = self.get_radius()
232 if radius is None:
233 lat = random_lat(rstate)
234 lon = rstate.uniform(-180., 180.)
235 elif radius == 0.0:
236 lat, lon = self.get_center_latlon()
237 else:
238 lat_center, lon_center = self.get_center_latlon()
239 while True:
240 north = rstate.uniform(-radius, radius)
241 east = rstate.uniform(-radius, radius)
242 if math.sqrt(north**2 + east**2) <= radius:
243 break
245 lat, lon = od.ne_to_latlon(lat_center, lon_center, north, east)
247 if not self.avoid_water or is_on_land(lat, lon):
248 logger.debug('location ok: %g, %g' % (lat, lon))
249 return lat, lon
251 if self.avoid_water:
252 sadd = ' (avoiding water)'
254 raise LocationGenerationError('Could not generate location%s.' % sadd)
256 def get_cartesian_offset(self, i):
257 rstate = self.get_rstate((i+1)*3+1)
258 north_shift = random_uniform(
259 rstate, self.north_shift_min, self.north_shift_max, 0.0)
260 east_shift = random_uniform(
261 rstate, self.east_shift_min, self.east_shift_max, 0.0)
263 return north_shift, east_shift
265 def get_depth(self, i):
266 rstate = self.get_rstate((i+1)*3+1)
267 return random_uniform(
268 rstate, self.depth_min, self.depth_max, 0.0)
270 def get_coordinates(self, i):
271 lat, lon = self.get_latlon(i)
272 north_shift, east_shift = self.get_cartesian_offset(i)
273 depth = self.get_depth(i)
274 return lat, lon, north_shift, east_shift, depth