Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/scenario/base.py: 72%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Base classes for the scenario generators. 

8''' 

9 

10import math 

11import hashlib 

12import logging 

13 

14import numpy as num 

15 

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 

20 

21logger = logging.getLogger('pyrocko.scenario.base') 

22 

23guts_prefix = 'pf.scenario' 

24 

25km = 1000. 

26d2r = num.pi/180. 

27N = 10000000 

28 

29coastlines = None 

30 

31 

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 

38 

39 

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. 

47 

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')) 

52 

53 return is_land 

54 

55 

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 

60 

61 

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 

69 

70 if avoid_water: 

71 sadd = ' (avoiding water)' 

72 

73 raise LocationGenerationError('Could not generate location%s.' % sadd) 

74 

75 

76def random_uniform(rstate, xmin, xmax, xdefault): 

77 if None in (xmin, xmax): 

78 return xdefault 

79 else: 

80 return rstate.uniform(xmin, xmax) 

81 

82 

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.') 

90 

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 

97 

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() 

104 

105 def clear(self): 

106 self._seed = None 

107 

108 def hash(self): 

109 return hashlib.sha1( 

110 (self.dump() + '\n\n%i' % self._retry_offset).encode('utf8'))\ 

111 .hexdigest() 

112 

113 def get_seed_offset(self): 

114 return int(self.hash(), base=16) % N 

115 

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) 

125 

126 def get_root(self): 

127 if self._parent is None: 

128 return self 

129 else: 

130 return self._parent.get_root() 

131 

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 

143 

144 return self._seed + self.get_seed_offset() 

145 

146 def get_rstate(self, i): 

147 return num.random.RandomState(int(self.get_seed() + i)) 

148 

149 def get_center_latlon(self): 

150 return self._parent.get_center_latlon() 

151 

152 def get_radius(self): 

153 return self._parent.get_radius() 

154 

155 def get_stations(self): 

156 return [] 

157 

158 

159class LocationGenerator(Generator): 

160 ''' 

161 Base class for generators providing random locations. 

162 ''' 

163 

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].') 

198 

199 def __init__(self, **kwargs): 

200 Generator.__init__(self, **kwargs) 

201 self._center_latlon = None 

202 

203 def clear(self): 

204 Generator.clear(self) 

205 self._center_latlon = None 

206 

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__) 

212 

213 if self._center_latlon is None: 

214 

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 

217 

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__) 

226 

227 return self._center_latlon 

228 

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 

236 

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 

254 

255 lat, lon = od.ne_to_latlon(lat_center, lon_center, north, east) 

256 

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 

260 

261 if self.avoid_water: 

262 sadd = ' (avoiding water)' 

263 

264 raise LocationGenerationError('Could not generate location%s.' % sadd) 

265 

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) 

272 

273 return north_shift, east_shift 

274 

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) 

279 

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