Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/dataset/crust2x2.py: 81%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5''' 

6Interface to use CRUST2.0 model by Laske, Masters and Reif. 

7 

8All functions defined in this module return SI units (m, m/s, kg/m^3). 

9 

10.. rubric:: Citation 

11 

12If you use this model, please refer to the 

13`REM web site <http://igppweb.ucsd.edu/~gabi/rem.html>`_ or cite: 

14 

15Bassin, C., Laske, G. and Masters, G., The Current Limits of Resolution for 

16Surface Wave Tomography in North America, EOS Trans AGU, 81, F897, 2000. A 

17description of CRUST 5.1 can be found in: Mooney, Laske and Masters, Crust 

185.1: a global crustal model at 5x5 degrees, JGR, 103, 727-747, 1998. 

19 

20.. rubric:: Usage 

21 

22:: 

23 

24 >>> from pyrocko import crust2x2 

25 >>> p = crust2x2.get_profile(10., 20.) 

26 >>> print p 

27 type, name: G2, Archean 0.5 km seds. 

28 elevation: 529 

29 crustal thickness: 38500 

30 average vp, vs, rho: 6460.7 3665.1 2867.5 

31 mantle ave. vp, vs, rho: 8200 4700 3400 

32 

33 0 3810 1940 920 ice 

34 0 1500 0 1020 water 

35 500 2500 1200 2100 soft sed. 

36 0 4000 2100 2400 hard sed. 

37 12500 6200 3600 2800 upper crust 

38 13000 6400 3600 2850 middle crust 

39 13000 6800 3800 2950 lower crust 

40 >>> print p.get_weeded() 

41 [[ 0. 500. 500. 13000. 13000. 26000. 26000. 39000. 39000.] 

42 [ 2500. 2500. 6200. 6200. 6400. 6400. 6800. 6800. 8200.] 

43 [ 1200. 1200. 3600. 3600. 3600. 3600. 3800. 3800. 4700.] 

44 [ 2100. 2100. 2800. 2800. 2850. 2850. 2950. 2950. 3400.]] 

45 

46 

47.. rubric:: Constants 

48 

49============== ============== 

50Layer id Layer name 

51============== ============== 

52LICE ice 

53LWATER water 

54LSOFTSED soft sediments 

55LHARDSED hard sediments 

56LUPPERCRUST upper crust 

57LMIDDLECRUST middle crust 

58LLOWERCRUST lower crust 

59LBELOWCRUST below crust 

60============== ============== 

61 

62''' 

63 

64import os 

65import copy 

66import math 

67from io import StringIO 

68 

69import numpy as num 

70 

71 

72LICE, LWATER, LSOFTSED, LHARDSED, LUPPERCRUST, LMIDDLECRUST, \ 

73 LLOWERCRUST, LBELOWCRUST = list(range(8)) 

74 

75 

76class Crust2Profile(object): 

77 ''' 

78 Representation of a CRUST2.0 key profile. 

79 ''' 

80 

81 layer_names = ( 

82 'ice', 'water', 'soft sed.', 'hard sed.', 'upper crust', 

83 'middle crust', 'lower crust', 'mantle') 

84 

85 def __init__(self, ident, name, vp, vs, rho, thickness, elevation): 

86 self._ident = ident 

87 self._name = name 

88 self._vp = vp 

89 self._vs = vs 

90 self._rho = rho 

91 self._thickness = thickness 

92 self._crustal_thickness = None 

93 self._elevation = elevation 

94 

95 def get_weeded(self, include_waterlayer=False): 

96 ''' 

97 Get layers used in the profile. 

98 

99 :param include_waterlayer: include water layer if ``True``. Default is 

100 ``False`` 

101 

102 :returns: NumPy array with rows ``depth``, ``vp``, ``vs``, ``density`` 

103 ''' 

104 depth = 0. 

105 layers = [] 

106 for ilayer, thickness, vp, vs, rho in zip( 

107 range(8), 

108 self._thickness, 

109 self._vp[:-1], 

110 self._vs[:-1], 

111 self._rho[:-1]): 

112 

113 if thickness == 0.0: 

114 continue 

115 

116 if not include_waterlayer and ilayer == LWATER: 

117 continue 

118 

119 layers.append([depth, vp, vs, rho]) 

120 layers.append([depth+thickness, vp, vs, rho]) 

121 depth += thickness 

122 

123 layers.append([ 

124 depth, 

125 self._vp[LBELOWCRUST], 

126 self._vs[LBELOWCRUST], 

127 self._rho[LBELOWCRUST]]) 

128 

129 return num.array(layers).T 

130 

131 def get_layer(self, ilayer): 

132 ''' 

133 Get parameters for a layer. 

134 

135 :param ilayer: id of layer 

136 :returns: thickness, vp, vs, density 

137 ''' 

138 

139 if ilayer == LBELOWCRUST: 

140 thickness = num.Inf 

141 else: 

142 thickness = self._thickness[ilayer] 

143 

144 return thickness, self._vp[ilayer], self._vs[ilayer], self._rho[ilayer] 

145 

146 def set_elevation(self, elevation): 

147 self._elevation = elevation 

148 

149 def set_layer_thickness(self, ilayer, thickness): 

150 self._thickness[ilayer] = thickness 

151 

152 def elevation(self): 

153 return self._elevation 

154 

155 def __str__(self): 

156 

157 vvp, vvs, vrho, vthi = self.averages() 

158 

159 return '''type, name: %s, %s 

160elevation: %15.5g 

161crustal thickness: %15.5g 

162average vp, vs, rho: %15.5g %15.5g %15.5g 

163mantle ave. vp, vs, rho: %15.5g %15.5g %15.5g 

164 

165%s''' % (self._ident, self._name, self._elevation, vthi, vvp, vvs, vrho, 

166 self._vp[LBELOWCRUST], self._vs[LBELOWCRUST], self._rho[LBELOWCRUST], 

167 '\n'.join([ 

168 '%15.5g %15.5g %15.5g %15.5g %s' % x 

169 for x in zip( 

170 self._thickness, self._vp[:-1], self._vs[:-1], self._rho[:-1], 

171 Crust2Profile.layer_names)])) 

172 

173 def crustal_thickness(self): 

174 ''' 

175 Get total crustal thickness. 

176 

177 Takes into account ice layer. 

178 Does not take into account water layer. 

179 ''' 

180 

181 return num.sum(self._thickness[3:]) + self._thickness[LICE] 

182 

183 def averages(self): 

184 ''' 

185 Get crustal averages for vp, vs, density and total crustal thickness. 

186 

187 Takes into account ice layer. 

188 Does not take into account water layer. 

189 ''' 

190 

191 vthi = self.crustal_thickness() 

192 vvp = num.sum(self._thickness[3:] / self._vp[3:-1]) + \ 

193 self._thickness[LICE] / self._vp[LICE] 

194 vvs = num.sum(self._thickness[3:] / self._vs[3:-1]) + \ 

195 self._thickness[LICE] / self._vs[LICE] 

196 vrho = num.sum(self._thickness[3:] * self._rho[3:-1]) + \ 

197 self._thickness[LICE] * self._rho[LICE] 

198 

199 vvp = vthi / vvp 

200 vvs = vthi / vvs 

201 vrho = vrho / vthi 

202 

203 return vvp, vvs, vrho, vthi 

204 

205 

206def _sa2arr(sa): 

207 return num.array([float(x) for x in sa], dtype=float) 

208 

209 

210def _wrap(x, mi, ma): 

211 if mi <= x and x <= ma: 

212 return x 

213 

214 return x - math.floor((x-mi)/(ma-mi)) * (ma-mi) 

215 

216 

217def _clip(x, mi, ma): 

218 return min(max(mi, x), ma) 

219 

220 

221class Crust2(object): 

222 ''' 

223 Access CRUST2.0 model. 

224 

225 :param directory: Directory with the data files which contain the 

226 CRUST2.0 model data. If this is set to ``None``, builtin CRUST2.0 

227 files are used. 

228 ''' 

229 

230 fn_keys = 'CNtype2_key.txt' 

231 fn_elevation = 'CNelevatio2.txt' 

232 fn_map = 'CNtype2.txt' 

233 

234 nlo = 180 

235 nla = 90 

236 

237 _instance = None 

238 

239 def __init__(self, directory=None): 

240 

241 self.profile_keys = [] 

242 self._directory = directory 

243 self._typemap = None 

244 self._load_crustal_model() 

245 

246 def get_profile(self, *args, **kwargs): 

247 ''' 

248 Get crustal profile at a specific location or raw profile for given 

249 key. 

250 

251 Get profile for location ``(lat, lon)``, or raw profile for given 

252 string key. 

253 

254 :rtype: instance of :py:class:`Crust2Profile` 

255 ''' 

256 

257 lat = kwargs.pop('lat', None) 

258 lon = kwargs.pop('lon', None) 

259 

260 if len(args) == 2: 

261 lat, lon = args 

262 

263 if lat is not None and lon is not None: 

264 return self._typemap[self._indices(float(lat), float(lon))] 

265 else: 

266 return self._raw_profiles[args[0]] 

267 

268 def _indices(self, lat, lon): 

269 lat = _clip(lat, -90., 90.) 

270 lon = _wrap(lon, -180., 180.) 

271 dlo = 360./Crust2.nlo 

272 dla = 180./Crust2.nla 

273 cola = 90.-lat 

274 ilat = _clip(int(cola/dla), 0, Crust2.nla-1) 

275 ilon = int((lon+180.)/dlo) % Crust2.nlo 

276 return ilat, ilon 

277 

278 def _load_crustal_model(self): 

279 

280 if self._directory is not None: 

281 path_keys = os.path.join(self._directory, Crust2.fn_keys) 

282 f = open(path_keys, 'r') 

283 else: 

284 from .crust2x2_data import decode, type2_key, type2, elevation 

285 

286 f = StringIO(decode(type2_key)) 

287 

288 # skip header 

289 for i in range(5): 

290 f.readline() 

291 

292 profiles = {} 

293 while True: 

294 line = f.readline() 

295 if not line: 

296 break 

297 ident, name = line.split(None, 1) 

298 line = f.readline() 

299 vp = _sa2arr(line.split()) * 1000. 

300 line = f.readline() 

301 vs = _sa2arr(line.split()) * 1000. 

302 line = f.readline() 

303 rho = _sa2arr(line.split()) * 1000. 

304 line = f.readline() 

305 toks = line.split() 

306 thickness = _sa2arr(toks[:-2]) * 1000. 

307 

308 assert ident not in profiles 

309 

310 profiles[ident] = Crust2Profile( 

311 ident.strip(), name.strip(), vp, vs, rho, thickness, 0.0) 

312 

313 f.close() 

314 

315 self._raw_profiles = profiles 

316 self.profile_keys = sorted(profiles.keys()) 

317 

318 if self._directory is not None: 

319 path_map = os.path.join(self._directory, Crust2.fn_map) 

320 f = open(path_map, 'r') 

321 else: 

322 f = StringIO(decode(type2)) 

323 

324 f.readline() # header 

325 

326 amap = {} 

327 for ila, line in enumerate(f): 

328 keys = line.split()[1:] 

329 for ilo, key in enumerate(keys): 

330 amap[ila, ilo] = copy.deepcopy(profiles[key]) 

331 

332 f.close() 

333 

334 if self._directory is not None: 

335 path_elevation = os.path.join(self._directory, Crust2.fn_elevation) 

336 f = open(path_elevation, 'r') 

337 

338 else: 

339 f = StringIO(decode(elevation)) 

340 

341 f.readline() 

342 for ila, line in enumerate(f): 

343 for ilo, s in enumerate(line.split()[1:]): 

344 p = amap[ila, ilo] 

345 p.set_elevation(float(s)) 

346 if p.elevation() < 0.: 

347 p.set_layer_thickness(LWATER, -p.elevation()) 

348 

349 f.close() 

350 

351 self._typemap = amap 

352 

353 @staticmethod 

354 def instance(): 

355 ''' 

356 Get the global default Crust2 instance. 

357 ''' 

358 

359 if Crust2._instance is None: 

360 Crust2._instance = Crust2() 

361 

362 return Crust2._instance 

363 

364 

365def get_profile_keys(): 

366 ''' 

367 Get list of all profile keys. 

368 ''' 

369 

370 crust2 = Crust2.instance() 

371 return list(crust2.profile_keys) 

372 

373 

374def get_profile(*args, **kwargs): 

375 ''' 

376 Get Crust2x2 profile for given location or profile key. 

377 

378 Get profile for (lat,lon) or raw profile for given string key. 

379 ''' 

380 

381 crust2 = Crust2.instance() 

382 return crust2.get_profile(*args, **kwargs) 

383 

384 

385def plot_crustal_thickness(crust2=None, filename='crustal_thickness.pdf'): 

386 ''' 

387 Create a quick and dirty plot of the crustal thicknesses defined in 

388 CRUST2.0. 

389 ''' 

390 

391 if crust2 is None: 

392 crust2 = Crust2.instance() 

393 

394 def func(lat, lon): 

395 return crust2.get_profile(lat, lon).crustal_thickness(), 

396 

397 plot(func, filename, zscaled_unit='km', zscaled_unit_factor=0.001) 

398 

399 

400def plot_vp_belowcrust(crust2=None, filename='vp_below_crust.pdf'): 

401 ''' 

402 Create a quick and dirty plot of vp below the crust, as defined in 

403 CRUST2.0. 

404 ''' 

405 

406 if crust2 is None: 

407 crust2 = Crust2.instance() 

408 

409 def func(lat, lon): 

410 return crust2.get_profile(lat, lon).get_layer(LBELOWCRUST)[1] 

411 

412 plot(func, filename, zscaled_unit='km/s', zscaled_unit_factor=0.001) 

413 

414 

415def plot(func, filename, **kwargs): 

416 nlats, nlons = 91, 181 

417 lats = num.linspace(-90., 90., nlats) 

418 lons = num.linspace(-180., 180., nlons) 

419 

420 vecfunc = num.vectorize(func, [float]) 

421 latss, lonss = num.meshgrid(lats, lons) 

422 thickness = vecfunc(latss, lonss) 

423 

424 from pyrocko.plot import gmtpy 

425 cm = gmtpy.cm 

426 marg = (1.5*cm, 2.5*cm, 1.5*cm, 1.5*cm) 

427 p = gmtpy.Simple( 

428 width=20*cm, height=10*cm, margins=marg, 

429 with_palette=True, **kwargs) 

430 

431 p.density_plot(gmtpy.tabledata(lons, lats, thickness.T)) 

432 p.save(filename)