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

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`