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

440 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 15:01 +0000

1# https://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5''' 

6Access to the `USGS Global Crustal Database 

7<https://earthquake.usgs.gov/data/crust>`_. 

8 

9The USGS Global Crustal Database provides empirical velocity measurements of 

10the earth for statistical analysis. 

11 

12.. rubric:: Citation 

13 

14Mooney, W. D., Laske, G., & Masters, T. G. (1998). CRUST 5.1: A global crustal 

15model at 5x5. Journal of Geophysical Research: Solid Earth, 103(B1), 727-747. 

16 

17''' 

18 

19import numpy as num 

20import copy 

21import logging 

22from os import path 

23 

24from pyrocko.guts import Object, String, Float, Int 

25from pyrocko.guts_array import Array 

26 

27from pyrocko.cake import LayeredModel, Material 

28from pyrocko.plot.cake_plot import my_model_plot, xscaled, yscaled 

29 

30from .crustdb_abbr import ageKey, provinceKey, referenceKey, pubYear # noqa 

31 

32logger = logging.getLogger('pyrocko.dataset.crustdb') 

33THICKNESS_HALFSPACE = 2 

34 

35db_url = 'https://mirror.pyrocko.org/gsc20130501.txt' 

36km = 1e3 

37vel_labels = { 

38 'vp': '$V_P$', 

39 'p': '$V_P$', 

40 'vs': '$V_S$', 

41 's': '$V_S$', 

42} 

43 

44 

45class DatabaseError(Exception): 

46 pass 

47 

48 

49class ProfileEmpty(Exception): 

50 pass 

51 

52 

53def _getCanvas(axes): 

54 

55 import matplotlib.pyplot as plt 

56 

57 if axes is None: 

58 fig = plt.figure() 

59 return fig, fig.gca() 

60 return axes.figure, axes 

61 

62 

63def xoffset_scale(offset, scale, ax): 

64 from matplotlib.ticker import ScalarFormatter, AutoLocator 

65 

66 class FormatVelocities(ScalarFormatter): 

67 @staticmethod 

68 def __call__(value, pos): 

69 return u'%.1f' % ((value-offset) * scale) 

70 

71 class OffsetLocator(AutoLocator): 

72 def tick_values(self, vmin, vmax): 

73 return [v + offset for v in 

74 AutoLocator.tick_values(self, vmin, vmax)] 

75 

76 ax.get_xaxis().set_major_formatter(FormatVelocities()) 

77 ax.get_xaxis().set_major_locator(OffsetLocator()) 

78 

79 

80class VelocityProfile(Object): 

81 uid = Int.T( 

82 optional=True, 

83 help='Unique ID of measurement') 

84 

85 lat = Float.T( 

86 help='Latitude [deg]') 

87 lon = Float.T( 

88 help='Longitude [deg]') 

89 elevation = Float.T( 

90 default=num.nan, 

91 help='Elevation [m]') 

92 vp = Array.T( 

93 shape=(None, 1), 

94 help='P Wave velocities [m/s]') 

95 vs = Array.T( 

96 shape=(None, 1), 

97 help='S Wave velocities [m/s]') 

98 d = Array.T( 

99 shape=(None, 1), 

100 help='Interface depth, top [m]') 

101 h = Array.T( 

102 shape=(None, 1), 

103 help='Interface thickness [m]') 

104 

105 heatflow = Float.T( 

106 optional=True, 

107 help='Heatflow [W/m^2]') 

108 geographical_location = String.T( 

109 optional=True, 

110 help='Geographic Location') 

111 geological_province = String.T( 

112 optional=True, 

113 help='Geological Province') 

114 geological_age = String.T( 

115 optional=True, 

116 help='Geological Age') 

117 measurement_method = Int.T( 

118 optional=True, 

119 help='Measurement method') 

120 publication_reference = String.T( 

121 optional=True, 

122 help='Publication Reference') 

123 publication_year__ = Int.T( 

124 help='Publication Date') 

125 

126 def __init__(self, *args, **kwargs): 

127 Object.__init__(self, *args, **kwargs) 

128 

129 self.h = num.abs(self.d - num.roll(self.d, -1)) 

130 self.h[-1] = 0 

131 self.nlayers = self.h.size 

132 

133 self.geographical_location = '%s (%s)' % ( 

134 provinceKey(self.geographical_location), 

135 self.geographical_location) 

136 

137 self.vs[self.vs == 0] = num.nan 

138 self.vp[self.vp == 0] = num.nan 

139 

140 self._step_vp = num.repeat(self.vp, 2) 

141 self._step_vs = num.repeat(self.vs, 2) 

142 self._step_d = num.roll(num.repeat(self.d, 2), -1) 

143 self._step_d[-1] = self._step_d[-2] + THICKNESS_HALFSPACE 

144 

145 @property 

146 def publication_year__(self): 

147 return pubYear(self.publication_reference) 

148 

149 def interpolateProfile(self, depths, phase='p', stepped=True): 

150 ''' 

151 Get a continuous velocity function at arbitrary depth. 

152 

153 :param depth: Depths to interpolate 

154 :type depth: :class:`numpy.ndarray` 

155 :param phase: P or S wave velocity, **p** or **s** 

156 :type phase: str 

157 :param stepped: Use a stepped velocity function or gradient 

158 :type stepped: bool 

159 :returns: velocities at requested depths 

160 :rtype: :class:`numpy.ndarray` 

161 ''' 

162 

163 if phase not in ['s', 'p']: 

164 raise AttributeError("Phase has to be either 'p' or 's'.") 

165 

166 if phase == 'p': 

167 vel = self._step_vp if stepped else self.vp 

168 elif phase == 's': 

169 vel = self._step_vs if stepped else self.vs 

170 d = self._step_d if stepped else self.d 

171 

172 if vel.size == 0: 

173 raise ProfileEmpty('Phase %s does not contain velocities' % phase) 

174 

175 try: 

176 res = num.interp(depths, d, vel, 

177 left=num.nan, right=num.nan) 

178 except ValueError: 

179 raise ValueError('Could not interpolate velocity profile.') 

180 

181 return res 

182 

183 def plot(self, axes=None): 

184 ''' 

185 Plot the velocity profile, see :class:`pyrocko.cake`. 

186 

187 :param axes: Axes to plot into. 

188 :type axes: :class:`matplotlib.axes.Axes` 

189 ''' 

190 

191 import matplotlib.pyplot as plt 

192 

193 fig, ax = _getCanvas(axes) 

194 my_model_plot(self.getLayeredModel(), axes=axes) 

195 ax.set_title('Global Crustal Database\n' 

196 'Velocity Structure at {p.lat:.4f}N, ' 

197 ' {p.lat:.4f}E (uid {p.uid})'.format(p=self)) 

198 if axes is None: 

199 plt.show() 

200 

201 def getLayeredModel(self): 

202 ''' 

203 Get a layered model, see :class:`pyrocko.cake.LayeredModel`. 

204 ''' 

205 def iterLines(): 

206 for il, m in enumerate(self.iterLayers()): 

207 yield self.d[il], m, '' 

208 

209 return LayeredModel.from_scanlines(iterLines()) 

210 

211 def iterLayers(self): 

212 ''' 

213 Iterator yielding a :class:`pyrocko.cake.Material` for each layer. 

214 ''' 

215 for il in range(self.nlayers): 

216 yield Material(vp=self.vp[il], 

217 vs=self.vs[il]) 

218 

219 @property 

220 def geog_loc_long(self): 

221 return provinceKey(self.geog_loc) 

222 

223 @property 

224 def geol_age_long(self): 

225 return ageKey(self.geol_age) 

226 

227 @property 

228 def has_s(self): 

229 return num.any(self.vp) 

230 

231 @property 

232 def has_p(self): 

233 return num.any(self.vs) 

234 

235 def _csv(self): 

236 output = '' 

237 for d in range(len(self.h)): 

238 output += ( 

239 '{p.uid}, {p.lat}, {p.lon},' 

240 ' {vp}, {vs}, {h}, {d}, {p.publication_reference}\n' 

241 ).format( 

242 p=self, 

243 vp=self.vp[d], vs=self.vs[d], h=self.h[d], d=self.d[d]) 

244 return output 

245 

246 

247class CrustDB(object): 

248 ''' 

249 CrustDB is a container for :class:`VelocityProfile` and provides 

250 functions for spatial selection, querying, processing and visualising 

251 data from the Global Crustal Database. 

252 ''' 

253 

254 def __init__(self, database_file=None, parent=None): 

255 self.profiles = [] 

256 self._velocity_matrix_cache = {} 

257 self.data_matrix = None 

258 self.name = None 

259 self.database_file = database_file 

260 

261 if parent is not None: 

262 pass 

263 elif database_file is not None: 

264 self._read(database_file) 

265 else: 

266 self._read(self._getRepositoryDatabase()) 

267 

268 def __len__(self): 

269 return len(self.profiles) 

270 

271 def __setitem__(self, key, value): 

272 if not isinstance(value, VelocityProfile): 

273 raise TypeError('Element is not a VelocityProfile') 

274 self.profiles[key] = value 

275 

276 def __delitem__(self, key): 

277 self.profiles.remove(key) 

278 

279 def __getitem__(self, key): 

280 return self.profiles[key] 

281 

282 def __str__(self): 

283 rstr = 'Container contains %d velocity profiles:\n\n' % self.nprofiles 

284 return rstr 

285 

286 @property 

287 def nprofiles(self): 

288 return len(self.profiles) 

289 

290 def append(self, value): 

291 if not isinstance(value, VelocityProfile): 

292 raise TypeError('Element is not a VelocityProfile') 

293 self.profiles.append(value) 

294 

295 def copy(self): 

296 return copy.deepcopy(self) 

297 

298 def lats(self): 

299 return num.array( 

300 [p.lat for p in self.profiles]) 

301 

302 def lons(self): 

303 return num.array( 

304 [p.lon for p in self.profiles]) 

305 

306 def _dataMatrix(self): 

307 if self.data_matrix is not None: 

308 return self.data_matrix 

309 

310 self.data_matrix = num.core.records.fromarrays( 

311 num.vstack([ 

312 num.concatenate([p.vp for p in self.profiles]), 

313 num.concatenate([p.vs for p in self.profiles]), 

314 num.concatenate([p.h for p in self.profiles]), 

315 num.concatenate([p.d for p in self.profiles]) 

316 ]), 

317 names='vp, vs, h, d') 

318 return self.data_matrix 

319 

320 def velocityMatrix(self, depth_range=(0, 60000.), ddepth=100., phase='p'): 

321 ''' 

322 Create a regular sampled velocity matrix. 

323 

324 :param depth_range: Depth range, ``(dmin, dmax)``, 

325 defaults to ``(0, 6000.)`` 

326 :type depth_range: tuple 

327 :param ddepth: Stepping in [m], defaults to ``100.`` 

328 :type ddepth: float 

329 :param phase: Phase to calculate ``p`` or ``s``, 

330 defaults to ``p`` 

331 :type phase: str 

332 :returns: Sample depths, veloctiy matrix 

333 :rtype: tuple, (sample_depth, :class:`numpy.ndarray`) 

334 ''' 

335 dmin, dmax = depth_range 

336 uid = '.'.join(map(repr, (dmin, dmax, ddepth, phase))) 

337 sdepth = num.linspace(dmin, dmax, int((dmax - dmin) / ddepth)) 

338 ndepth = sdepth.size 

339 

340 if uid not in self._velocity_matrix_cache: 

341 vel_mat = num.empty((self.nprofiles, ndepth)) 

342 for ip, profile in enumerate(self.profiles): 

343 vel_mat[ip, :] = profile.interpolateProfile(sdepth, 

344 phase=phase) 

345 self._velocity_matrix_cache[uid] = num.ma.masked_invalid(vel_mat) 

346 

347 return sdepth, self._velocity_matrix_cache[uid] 

348 

349 def rmsRank(self, ref_profile, depth_range=(0, 3500.), ddepth=100., 

350 phase='p'): 

351 ''' 

352 Correlates ``ref_profile`` to each profile in the database. 

353 

354 :param ref_profile: Reference profile 

355 :type ref_profile: :class:`VelocityProfile` 

356 :param depth_range: Depth range in [m], ``(dmin, dmax)``, 

357 defaults to ``(0, 35000.)`` 

358 :type depth_range: tuple 

359 :param ddepth: Stepping in [m], defaults to ``100.`` 

360 :type ddepth: float 

361 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p`` 

362 :type phase: str 

363 :returns: RMS factor length of N_profiles 

364 :rtype: :class:`numpy.ndarray` 

365 ''' 

366 if not isinstance(ref_profile, VelocityProfile): 

367 raise ValueError('ref_profile is not a VelocityProfile') 

368 

369 sdepth, vel_matrix = self.velocityMatrix(depth_range, ddepth, 

370 phase=phase) 

371 ref_vel = ref_profile.interpolateProfile(sdepth, phase=phase) 

372 

373 rms = num.empty(self.nprofiles) 

374 for p in range(self.nprofiles): 

375 profile = vel_matrix[p, :] 

376 rms[p] = num.sqrt(profile**2 - ref_vel**2).sum() / ref_vel.size 

377 return rms 

378 

379 def histogram2d(self, depth_range=(0., 60000.), vel_range=None, 

380 ddepth=100., dvbin=100., ddbin=2000., phase='p'): 

381 ''' 

382 Create a 2D Histogram of all the velocity profiles. 

383 

384 Check :func:`numpy.histogram2d` for more information. 

385 

386 :param depth_range: Depth range in [m], ``(dmin, dmax)``, 

387 defaults to ``(0., 60000.)`` 

388 :type depth_range: tuple 

389 :param vel_range: Depth range, ``(vmin, vmax)``, 

390 defaults to ``(5500., 8500.)`` 

391 :type vel_range: tuple 

392 :param ddepth: Stepping in [km], defaults to ``100.`` 

393 :type ddepth: float 

394 :param dvbin: Bin size in velocity dimension [m/s], defaults to 100. 

395 :type dvbin: float 

396 :param dvbin: Bin size in depth dimension [m], defaults to 2000. 

397 :type dvbin: float 

398 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p`` 

399 :type phase: str 

400 

401 :returns: :func:`numpy.histogram2d` 

402 :rtype: tuple 

403 ''' 

404 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth, phase=phase) 

405 d_vec = num.tile(sdepth, self.nprofiles) 

406 

407 # Velocity and depth bins 

408 if vel_range is None: 

409 vel_range = ((v_mat.min() // 1e2) * 1e2, 

410 (v_mat.max() // 1e2) * 1e2) 

411 nvbins = int((vel_range[1] - vel_range[0]) / dvbin) 

412 ndbins = int((depth_range[1] - depth_range[0]) / ddbin) 

413 

414 return num.histogram2d(v_mat.flatten(), d_vec, 

415 range=(vel_range, depth_range), 

416 bins=[nvbins, ndbins], 

417 density=False) 

418 

419 def meanVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'): 

420 ''' 

421 Get mean and standard deviation of velocity profile. 

422 

423 :param depth_range: Depth range in [m], ``(dmin, dmax)``, 

424 defaults to ``(0., 60000.)`` 

425 :type depth_range: tuple 

426 :param ddepth: Stepping in [m], defaults to ``100.`` 

427 :type ddepth: float 

428 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p`` 

429 :type phase: str 

430 :returns: depth vector, mean velocities, standard deviations 

431 :rtype: :py:class:`tuple` of :class:`numpy.ndarray` 

432 ''' 

433 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth, phase=phase) 

434 v_mean = num.ma.mean(v_mat, axis=0) 

435 v_std = num.ma.std(v_mat, axis=0) 

436 

437 return sdepth, v_mean.flatten(), v_std.flatten() 

438 

439 def modeVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'): 

440 ''' 

441 Get mode velocity profile and standard deviation. 

442 

443 :param depth_range: Depth range in [m], ``(dmin, dmax)``, 

444 defaults to ``(0., 60000.)`` 

445 :type depth_range: tuple 

446 :param ddepth: Stepping in [m], defaults to ``100.`` 

447 :type ddepth: float 

448 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p`` 

449 :type phase: str 

450 :returns: depth vector, mode velocity, number of counts at each depth 

451 :rtype: :py:class:`tuple` of :class:`numpy.ndarray` 

452 ''' 

453 import scipy.stats 

454 

455 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth, phase=phase) 

456 v_mode, v_counts = scipy.stats.mstats.mode(v_mat, axis=0) 

457 return sdepth, v_mode.flatten(), v_counts.flatten() 

458 

459 def medianVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'): 

460 ''' 

461 Median velocity profile plus std variation. 

462 

463 :param depth_range: Depth range in [m], ``(dmin, dmax)``, 

464 defaults to ``(0., 60000.)`` 

465 :type depth_range: tuple 

466 :param ddepth: Stepping in [m], defaults to ``100.`` 

467 :type ddepth: float 

468 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p`` 

469 :type phase: str 

470 :returns: depth vector, median velocities, standard deviations 

471 :rtype: :py:class:`tuple` of :class:`numpy.ndarray` 

472 ''' 

473 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth, phase=phase) 

474 v_mean = num.ma.median(v_mat, axis=0) 

475 v_std = num.ma.std(v_mat, axis=0) 

476 

477 return sdepth, v_mean.flatten(), v_std.flatten() 

478 

479 def plotHistogram(self, vel_range=None, bins=36, phase='vp', 

480 axes=None): 

481 ''' 

482 Plot 1D histogram of seismic velocities in the container. 

483 

484 :param vel_range: Velocity range, defaults to (5.5, 8.5) 

485 :type vel_range: tuple 

486 :param bins: bins, defaults to 30 (see :func:`numpy.histogram`) 

487 :type bins: int 

488 :param phase: Property to plot out of ``['vp', 'vs']``, 

489 defaults to 'vp' 

490 :type phase: str 

491 :param figure: Figure to plot in, defaults to None 

492 :type figure: :class:`matplotlib.figure.Figure` 

493 ''' 

494 

495 import matplotlib.pyplot as plt 

496 

497 fig, ax = _getCanvas(axes) 

498 

499 if phase not in ['vp', 'vs']: 

500 raise AttributeError('phase has to be either vp or vs') 

501 

502 data = self._dataMatrix()[phase] 

503 

504 ax.hist(data, weights=self.data_matrix['h'], 

505 range=vel_range, bins=bins, 

506 color='g', alpha=.5) 

507 ax.text(.95, .95, '%d Profiles' % self.nprofiles, 

508 transform=ax.transAxes, fontsize=10, 

509 va='top', ha='right', alpha=.7) 

510 

511 ax.set_title('Distribution of %s' % vel_labels[phase]) 

512 ax.set_xlabel('%s [km/s]' % vel_labels[phase]) 

513 ax.set_ylabel('Cumulative occurrence [N]') 

514 xscaled(1./km, ax) 

515 ax.yaxis.grid(alpha=.4) 

516 

517 if self.name is not None: 

518 ax.set_title('%s for %s' % (ax.get_title(), self.name)) 

519 

520 if axes is None: 

521 plt.show() 

522 

523 def plot(self, depth_range=(0, 60000.), ddepth=100., ddbin=2000., 

524 vel_range=None, dvbin=100., 

525 percent=False, 

526 plot_mode=True, plot_median=True, plot_mean=False, 

527 show_cbar=True, 

528 aspect=.02, 

529 phase='p', 

530 axes=None): 

531 ''' 

532 Plot a two 2D Histogram of seismic velocities. 

533 

534 :param depth_range: Depth range, ``(dmin, dmax)``, 

535 defaults to ``(0, 60)`` 

536 :type depth_range: tuple 

537 :param vel_range: Velocity range, ``(vmin, vmax)`` 

538 :type vel_range: tuple 

539 :param ddepth: Stepping in [m], defaults to ``.1`` 

540 :type ddepth: float 

541 :param dvbin: Bin size in velocity dimension [m/s], defaults to .1 

542 :type dvbin: float 

543 :param dvbin: Bin size in depth dimension [m], defaults to 2000. 

544 :type dvbin: float 

545 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p`` 

546 :type phase: str 

547 :param plot_mode: Plot the Mode 

548 :type plot_mode: bool 

549 :param plot_mean: Plot the Mean 

550 :type plot_mean: bool 

551 :param plot_median: Plot the Median 

552 :type plot_median: bool 

553 :param axes: Axes to plot into, defaults to None 

554 :type axes: :class:`matplotlib.axes.Axes` 

555 ''' 

556 

557 import matplotlib.pyplot as plt 

558 

559 fig, ax = _getCanvas(axes) 

560 

561 ax = fig.gca() 

562 

563 if vel_range is not None: 

564 vmin, vmax = vel_range 

565 dmin, dmax = depth_range 

566 

567 vfield, vedg, dedg = self.histogram2d(vel_range=vel_range, 

568 depth_range=depth_range, 

569 ddepth=ddepth, dvbin=dvbin, 

570 ddbin=ddbin, phase=phase) 

571 vfield /= (ddbin / ddepth) 

572 

573 if percent: 

574 vfield /= vfield.sum(axis=1)[num.newaxis, :] 

575 

576 grid_ext = [vedg[0], vedg[-1], dedg[-1], dedg[0]] 

577 histogram = ax.imshow(vfield.swapaxes(0, 1), 

578 interpolation='nearest', 

579 extent=grid_ext, aspect=aspect) 

580 

581 if show_cbar: 

582 cticks = num.unique( 

583 num.arange(0, vfield.max(), vfield.max() // 10).round()) 

584 cbar = fig.colorbar(histogram, ticks=cticks, format='%1i', 

585 orientation='horizontal') 

586 if percent: 

587 cbar.set_label('Percent') 

588 else: 

589 cbar.set_label('Number of Profiles') 

590 

591 if plot_mode: 

592 sdepth, vel_mode, _ = self.modeVelocity(depth_range=depth_range, 

593 ddepth=ddepth, phase=phase) 

594 ax.plot(vel_mode[sdepth < dmax] + ddepth/2, 

595 sdepth[sdepth < dmax], 

596 alpha=.8, color='w', label='Mode') 

597 

598 if plot_mean: 

599 sdepth, vel_mean, _ = self.meanVelocity(depth_range=depth_range, 

600 ddepth=ddepth, phase=phase) 

601 ax.plot(vel_mean[sdepth < dmax] + ddepth/2, 

602 sdepth[sdepth < dmax], 

603 alpha=.8, color='w', linestyle='--', label='Mean') 

604 

605 if plot_median: 

606 sdepth, vel_median, _ = self.medianVelocity( 

607 depth_range=depth_range, 

608 ddepth=ddepth, phase=phase) 

609 ax.plot(vel_median[sdepth < dmax] + ddepth/2, 

610 sdepth[sdepth < dmax], 

611 alpha=.8, color='w', linestyle=':', label='Median') 

612 

613 ax.grid(True, which='both', color='w', linewidth=.8, alpha=.4) 

614 

615 ax.text(.025, .025, '%d Profiles' % self.nprofiles, 

616 color='w', alpha=.7, 

617 transform=ax.transAxes, fontsize=9, va='bottom', ha='left') 

618 

619 ax.set_title('Crustal Velocity Distribution') 

620 ax.set_xlabel('%s [km/s]' % vel_labels[phase]) 

621 ax.set_ylabel('Depth [km]') 

622 yscaled(1./km, ax) 

623 xoffset_scale(dvbin/2, 1./km, ax) 

624 ax.set_xlim(vel_range) 

625 

626 if self.name is not None: 

627 ax.set_title('%s for %s' % (ax.get_title(), self.name)) 

628 

629 if plot_mode or plot_mean or plot_median: 

630 leg = ax.legend(loc=1, fancybox=True, prop={'size': 10.}) 

631 leg.get_frame().set_alpha(.6) 

632 

633 if axes is None: 

634 plt.show() 

635 

636 def plotVelocitySurface(self, v_max, d_min=0., d_max=6000., axes=None): 

637 ''' 

638 Plot a triangulated a depth surface exceeding velocity. 

639 ''' 

640 

641 import matplotlib.pyplot as plt 

642 

643 fig, ax = _getCanvas(axes) 

644 d = self.exceedVelocity(v_max, d_min, d_max) 

645 lons = self.lons()[d > 0] 

646 lats = self.lats()[d > 0] 

647 d = d[d > 0] 

648 

649 ax.tricontourf(lats, lons, d) 

650 

651 if axes is None: 

652 plt.show() 

653 

654 def plotMap(self, outfile, **kwargs): 

655 from pyrocko.plot import gmtpy 

656 lats = self.lats() 

657 lons = self.lons() 

658 s, n, w, e = (lats.min(), lats.max(), lons.min(), lons.max()) 

659 

660 def darken(c, f=0.7): 

661 return (c[0]*f, c[1]*f, c[2]*f) 

662 

663 gmt = gmtpy.GMT() 

664 gmt.psbasemap(B='40/20', 

665 J='M0/12', 

666 R='%f/%f/%f/%f' % (w, e, s, n)) 

667 gmt.pscoast(R=True, J=True, 

668 D='i', S='216/242/254', A=10000, 

669 W='.2p') 

670 gmt.psxy(R=True, J=True, 

671 in_columns=[lons, lats], 

672 S='c2p', G='black') 

673 gmt.save(outfile) 

674 

675 def exceedVelocity(self, v_max, d_min=0, d_max=60): 

676 ''' 

677 Returns the last depth ``v_max`` has not been exceeded. 

678 

679 :param v_max: maximal velocity 

680 :type vmax: float 

681 :param dz: depth is sampled in dz steps 

682 :type dz: float 

683 :param d_max: maximum depth 

684 :type d_max: int 

685 :param d_min: minimum depth 

686 :type d_min: int 

687 :returns: Lat, Lon, Depth and uid where ``v_max`` is exceeded 

688 :rtype: :py:class:`list` of :py:class:`numpy.ndarray` 

689 ''' 

690 self.profile_exceed_velocity = num.empty(len(self.profiles)) 

691 self.profile_exceed_velocity[:] = num.nan 

692 

693 for ip, profile in enumerate(self.profiles): 

694 for il in range(len(profile.d)): 

695 if profile.d[il] <= d_min\ 

696 or profile.d[il] >= d_max: 

697 continue 

698 if profile.vp[il] < v_max: 

699 continue 

700 else: 

701 self.profile_exceed_velocity[ip] = profile.d[il] 

702 break 

703 return self.profile_exceed_velocity 

704 

705 def selectRegion(self, west, east, south, north): 

706 ''' 

707 Select profiles within a region by geographic corner coordinates. 

708 

709 :param west: west corner 

710 :type west: float 

711 :param east: east corner 

712 :type east: float 

713 :param south: south corner 

714 :type south: float 

715 :param north: north corner 

716 :type north: float 

717 :returns: Selected profiles 

718 :rtype: :class:`CrustDB` 

719 ''' 

720 r_container = self._emptyCopy() 

721 

722 for profile in self.profiles: 

723 if profile.lon >= west and profile.lon <= east \ 

724 and profile.lat <= north and profile.lat >= south: 

725 r_container.append(profile) 

726 

727 return r_container 

728 

729 def selectPolygon(self, poly): 

730 ''' 

731 Select profiles within a polygon. 

732 

733 The algorithm is called the **Ray Casting Method** 

734 

735 :param poly: Latitude Longitude pairs of the polygon 

736 :type param: list of :class:`numpy.ndarray` 

737 :returns: Selected profiles 

738 :rtype: :class:`CrustDB` 

739 ''' 

740 r_container = self._emptyCopy() 

741 

742 for profile in self.profiles: 

743 x = profile.lon 

744 y = profile.lat 

745 

746 inside = False 

747 p1x, p1y = poly[0] 

748 for p2x, p2y in poly: 

749 if y >= min(p1y, p2y): 

750 if y <= max(p1y, p2y): 

751 if x <= max(p1x, p2x): 

752 if p1y != p2y: 

753 xints = (y - p1y) * (p2x - p1x) / \ 

754 (p2y - p1y) + p1x 

755 if p1x == p2x or x <= xints: 

756 inside = not inside 

757 p1x, p1y = p2x, p2y 

758 if inside: 

759 r_container.append(profile) 

760 

761 return r_container 

762 

763 def selectLocation(self, lat, lon, radius=10): 

764 ''' 

765 Select profiles at a geographic location within a ``radius``. 

766 

767 :param lat: Latitude in [deg] 

768 :type lat: float 

769 :param lon: Longitude in [deg] 

770 :type lon: float 

771 :param radius: Radius in [deg] 

772 :type radius: float 

773 :returns: Selected profiles 

774 :rtype: :class:`CrustDB` 

775 ''' 

776 r_container = self._emptyCopy() 

777 logger.info('Selecting location %f, %f (r=%f)...' % (lat, lon, radius)) 

778 for profile in self.profiles: 

779 if num.sqrt((lat - profile.lat)**2 + 

780 (lon - profile.lon)**2) <= radius: 

781 r_container.append(profile) 

782 

783 return r_container 

784 

785 def selectMinLayers(self, nlayers): 

786 ''' 

787 Select profiles with more than ``nlayers``. 

788 

789 :param nlayers: Minimum number of layers 

790 :type nlayers: int 

791 :returns: Selected profiles 

792 :rtype: :class:`CrustDB` 

793 ''' 

794 r_container = self._emptyCopy() 

795 logger.info('Selecting minimum %d layers...' % nlayers) 

796 

797 for profile in self.profiles: 

798 if profile.nlayers >= nlayers: 

799 r_container.append(profile) 

800 

801 return r_container 

802 

803 def selectMaxLayers(self, nlayers): 

804 ''' 

805 Select profiles with more than ``nlayers``. 

806 

807 :param nlayers: Maximum number of layers 

808 :type nlayers: int 

809 :returns: Selected profiles 

810 :rtype: :class:`CrustDB` 

811 ''' 

812 r_container = self._emptyCopy() 

813 logger.info('Selecting maximum %d layers...' % nlayers) 

814 

815 for profile in self.profiles: 

816 if profile.nlayers <= nlayers: 

817 r_container.append(profile) 

818 

819 return r_container 

820 

821 def selectMinDepth(self, depth): 

822 ''' 

823 Select profiles describing layers deeper than ``depth``. 

824 

825 :param depth: Minumum depth in [m] 

826 :type depth: float 

827 :returns: Selected profiles 

828 :rtype: :class:`CrustDB` 

829 ''' 

830 r_container = self._emptyCopy() 

831 logger.info('Selecting minimum depth %f m...' % depth) 

832 

833 for profile in self.profiles: 

834 if profile.d.max() >= depth: 

835 r_container.append(profile) 

836 return r_container 

837 

838 def selectMaxDepth(self, depth): 

839 ''' 

840 Select profiles describing layers shallower than ``depth``. 

841 

842 :param depth: Maximum depth in [m] 

843 :type depth: float 

844 :returns: Selected profiles 

845 :rtype: :class:`CrustDB` 

846 ''' 

847 r_container = self._emptyCopy() 

848 logger.info('Selecting maximum depth %f m...' % depth) 

849 

850 for profile in self.profiles: 

851 if profile.d.max() <= depth: 

852 r_container.append(profile) 

853 return r_container 

854 

855 def selectVp(self): 

856 ''' 

857 Select profiles describing P wave velocity. 

858 

859 :returns Selected profiles 

860 :rtype: :class:`CrustDB` 

861 ''' 

862 r_container = self._emptyCopy() 

863 logger.info('Selecting profiles providing Vp...') 

864 

865 for profile in self.profiles: 

866 if not num.all(num.isnan(profile.vp)): 

867 r_container.append(profile) 

868 return r_container 

869 

870 def selectVs(self): 

871 ''' 

872 Select profiles describing P wave velocity. 

873 

874 :returns: Selected profiles 

875 :rtype: :class:`CrustDB` 

876 ''' 

877 r_container = self._emptyCopy() 

878 logger.info('Selecting profiles providing Vs...') 

879 

880 for profile in self.profiles: 

881 if not num.all(num.isnan(profile.vs)): 

882 r_container.append(profile) 

883 return r_container 

884 

885 def _emptyCopy(self): 

886 r_container = CrustDB(parent=self) 

887 r_container.name = self.name 

888 return r_container 

889 

890 def exportCSV(self, filename=None): 

891 ''' 

892 Export as CSV file. 

893 

894 :param filename: Export filename 

895 :type filename: str 

896 ''' 

897 with open(filename, 'w') as file: 

898 file.write('# uid, Lat, Lon, vp, vs, H, Depth, Reference\n') 

899 for profile in self.profiles: 

900 file.write(profile._csv()) 

901 

902 def exportYAML(self, filename=None): 

903 ''' 

904 Export as YAML file. 

905 

906 :param filename: Export filename 

907 :type filename: str 

908 ''' 

909 with open(filename, 'w') as file: 

910 for profile in self.profiles: 

911 file.write(profile.__str__()) 

912 

913 @classmethod 

914 def readDatabase(cls, database_file): 

915 db = cls() 

916 CrustDB._read(db, database_file) 

917 return db 

918 

919 @staticmethod 

920 def _getRepositoryDatabase(): 

921 from pyrocko import config 

922 

923 name = path.basename(db_url) 

924 data_path = path.join(config.config().crustdb_dir, name) 

925 if not path.exists(data_path): 

926 from pyrocko import util 

927 util.download_file(db_url, data_path, None, None) 

928 

929 return data_path 

930 

931 def _read(self, database_file): 

932 ''' 

933 Reads in the the GSN database file and puts it in CrustDB. 

934 

935 File format: 

936 

937 uid lat/lon vp vs hc depth 

938 2 29.76N 2.30 .00 2.00 .00 s 25.70 .10 .00 NAC-CO 5 U 

939 96.31W 3.94 .00 5.30 2.00 s 33.00 MCz 39.00 61C.3 EXC 

940 5.38 .00 12.50 7.30 c 

941 6.92 .00 13.20 19.80 c 

942 8.18 .00 .00 33.00 m 

943 

944 3 34.35N 3.00 .00 3.00 .00 s 35.00 1.60 .00 NAC-BR 4 R 

945 117.83W 6.30 .00 16.50 3.00 38.00 MCz 55.00 63R.1 ORO 

946 7.00 .00 18.50 19.50 

947 7.80 .00 .00 38.00 m 

948 

949 

950 :param database_file: path to database file, type string 

951 

952 ''' 

953 

954 def get_empty_record(): 

955 meta = { 

956 'uid': num.nan, 

957 'geographical_location': None, 

958 'geological_province': None, 

959 'geological_age': None, 

960 'elevation': num.nan, 

961 'heatflow': num.nan, 

962 'measurement_method': None, 

963 'publication_reference': None 

964 } 

965 # vp, vs, h, d, lat, lon, meta 

966 return (num.zeros(128, dtype=num.float32), 

967 num.zeros(128, dtype=num.float32), 

968 num.zeros(128, dtype=num.float32), 

969 num.zeros(128, dtype=num.float32), 

970 0., 0., meta) 

971 

972 nlayers = [] 

973 

974 def add_record(vp, vs, h, d, lat, lon, meta, nlayer): 

975 if nlayer == 0: 

976 return 

977 self.append(VelocityProfile( 

978 vp=vp[:nlayer] * km, 

979 vs=vs[:nlayer] * km, 

980 h=h[:nlayer] * km, 

981 d=d[:nlayer] * km, 

982 lat=lat, lon=lon, 

983 **meta)) 

984 nlayers.append(nlayer) 

985 

986 vp, vs, h, d, lat, lon, meta = get_empty_record() 

987 ilayer = 0 

988 with open(database_file, 'r') as database: 

989 for line, dbline in enumerate(database): 

990 if dbline.isspace(): 

991 if not len(d) == 0: 

992 add_record(vp, vs, h, d, lat, lon, meta, ilayer) 

993 if not len(vp) == len(h): 

994 raise DatabaseError( 

995 'Inconsistent database, check line %d!\n\tDebug: ' 

996 % line, lat, lon, vp, vs, h, d, meta) 

997 

998 vp, vs, h, d, lat, lon, meta = get_empty_record() 

999 ilayer = 0 

1000 else: 

1001 try: 

1002 if ilayer == 0: 

1003 lat = float(dbline[8:13]) 

1004 if dbline[13] == b'S': 

1005 lat = -lat 

1006 # Additional meta data 

1007 meta['uid'] = int(dbline[0:6]) 

1008 meta['elevation'] = float(dbline[52:57]) 

1009 meta['heatflow'] = float(dbline[58:64]) 

1010 if meta['heatflow'] == 0.: 

1011 meta['heatflow'] = None 

1012 meta['geographical_location'] =\ 

1013 dbline[66:72].strip() 

1014 meta['measurement_method'] = dbline[77] 

1015 if ilayer == 1: 

1016 lon = float(dbline[7:13]) 

1017 if dbline[13] == b'W': 

1018 lon = -lon 

1019 # Additional meta data 

1020 meta['geological_age'] = dbline[54:58].strip() 

1021 meta['publication_reference'] =\ 

1022 dbline[66:72].strip() 

1023 meta['geological_province'] = dbline[74:78].strip() 

1024 try: 

1025 vp[ilayer] = dbline[17:21] 

1026 vs[ilayer] = dbline[23:27] 

1027 h[ilayer] = dbline[28:34] 

1028 d[ilayer] = dbline[35:41] 

1029 except ValueError: 

1030 pass 

1031 except ValueError: 

1032 logger.warning( 

1033 'Could not interpret line %d, skipping\n%s' % 

1034 (line, dbline)) 

1035 while not database.readlines(): 

1036 pass 

1037 vp, vs, h, d, lat, lon, meta = get_empty_record() 

1038 ilayer += 1 

1039 # Append last profile 

1040 add_record(vp, vs, h, d, lat, lon, meta, ilayer) 

1041 logger.info('Loaded %d profiles from %s' % 

1042 (self.nprofiles, database_file))