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