1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5''' 

6Access to the USGS Global Crustal Database. 

7 

8Simple queries and statistical analysis 

9''' 

10 

11import numpy as num 

12import copy 

13import logging 

14from os import path 

15 

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

17from pyrocko.guts_array import Array 

18 

19from pyrocko.cake import LayeredModel, Material 

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

21 

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

23 

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

25THICKNESS_HALFSPACE = 2 

26 

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

28km = 1e3 

29vel_labels = { 

30 'vp': '$V_P$', 

31 'p': '$V_P$', 

32 'vs': '$V_S$', 

33 's': '$V_S$', 

34} 

35 

36 

37class DatabaseError(Exception): 

38 pass 

39 

40 

41class ProfileEmpty(Exception): 

42 pass 

43 

44 

45def _getCanvas(axes): 

46 

47 import matplotlib.pyplot as plt 

48 

49 if axes is None: 

50 fig = plt.figure() 

51 return fig, fig.gca() 

52 return axes.figure, axes 

53 

54 

55def xoffset_scale(offset, scale, ax): 

56 from matplotlib.ticker import ScalarFormatter, AutoLocator 

57 

58 class FormatVelocities(ScalarFormatter): 

59 @staticmethod 

60 def __call__(value, pos): 

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

62 

63 class OffsetLocator(AutoLocator): 

64 def tick_values(self, vmin, vmax): 

65 return [v + offset for v in 

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

67 

68 ax.get_xaxis().set_major_formatter(FormatVelocities()) 

69 ax.get_xaxis().set_major_locator(OffsetLocator()) 

70 

71 

72class VelocityProfile(Object): 

73 uid = Int.T( 

74 optional=True, 

75 help='Unique ID of measurement') 

76 

77 lat = Float.T( 

78 help='Latitude [deg]') 

79 lon = Float.T( 

80 help='Longitude [deg]') 

81 elevation = Float.T( 

82 default=num.nan, 

83 help='Elevation [m]') 

84 vp = Array.T( 

85 shape=(None, 1), 

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

87 vs = Array.T( 

88 shape=(None, 1), 

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

90 d = Array.T( 

91 shape=(None, 1), 

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

93 h = Array.T( 

94 shape=(None, 1), 

95 help='Interface thickness [m]') 

96 

97 heatflow = Float.T( 

98 optional=True, 

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

100 geographical_location = String.T( 

101 optional=True, 

102 help='Geographic Location') 

103 geological_province = String.T( 

104 optional=True, 

105 help='Geological Province') 

106 geological_age = String.T( 

107 optional=True, 

108 help='Geological Age') 

109 measurement_method = Int.T( 

110 optional=True, 

111 help='Measurement method') 

112 publication_reference = String.T( 

113 optional=True, 

114 help='Publication Reference') 

115 publication_year__ = Int.T( 

116 help='Publication Date') 

117 

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

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

120 

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

122 self.h[-1] = 0 

123 self.nlayers = self.h.size 

124 

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

126 provinceKey(self.geographical_location), 

127 self.geographical_location) 

128 

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

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

131 

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

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

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

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

136 

137 @property 

138 def publication_year__(self): 

139 return pubYear(self.publication_reference) 

140 

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

142 ''' 

143 Get a continuous velocity function at arbitrary depth. 

144 

145 :param depth: Depths to interpolate 

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

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

148 :type phase: str, optional 

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

150 :type stepped: bool 

151 :returns: velocities at requested depths 

152 :rtype: :class:`numpy.ndarray` 

153 ''' 

154 

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

156 raise AttributeError('Phase has to be either \'p\' or \'s\'.') 

157 

158 if phase == 'p': 

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

160 elif phase == 's': 

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

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

163 

164 if vel.size == 0: 

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

166 

167 try: 

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

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

170 except ValueError: 

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

172 

173 return res 

174 

175 def plot(self, axes=None): 

176 ''' 

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

178 

179 :param axes: Axes to plot into. 

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

181 ''' 

182 

183 import matplotlib.pyplot as plt 

184 

185 fig, ax = _getCanvas(axes) 

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

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

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

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

190 if axes is None: 

191 plt.show() 

192 

193 def getLayeredModel(self): 

194 ''' 

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

196 ''' 

197 def iterLines(): 

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

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

200 

201 return LayeredModel.from_scanlines(iterLines()) 

202 

203 def iterLayers(self): 

204 ''' 

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

206 ''' 

207 for il in range(self.nlayers): 

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

209 vs=self.vs[il]) 

210 

211 @property 

212 def geog_loc_long(self): 

213 return provinceKey(self.geog_loc) 

214 

215 @property 

216 def geol_age_long(self): 

217 return ageKey(self.geol_age) 

218 

219 @property 

220 def has_s(self): 

221 return num.any(self.vp) 

222 

223 @property 

224 def has_p(self): 

225 return num.any(self.vs) 

226 

227 def get_weeded(self): 

228 ''' 

229 Get weeded representation of layers used in the profile. 

230 See :func:`pyrocko.cake.get_weeded` for details. 

231 ''' 

232 weeded = num.zeros((self.nlayers, 4)) 

233 weeded[:, 0] = self.d 

234 weeded[:, 1] = self.vp 

235 weeded[:, 2] = self.vs 

236 

237 def _csv(self): 

238 output = '' 

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

240 output += ( 

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

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

243 ).format( 

244 p=self, 

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

246 return output 

247 

248 

249class CrustDB(object): 

250 ''' 

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

252 functions for spatial selection, querying, processing and visualising 

253 data from the Global Crustal Database. 

254 ''' 

255 

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

257 self.profiles = [] 

258 self._velocity_matrix_cache = {} 

259 self.data_matrix = None 

260 self.name = None 

261 self.database_file = database_file 

262 

263 if parent is not None: 

264 pass 

265 elif database_file is not None: 

266 self._read(database_file) 

267 else: 

268 self._read(self._getRepositoryDatabase()) 

269 

270 def __len__(self): 

271 return len(self.profiles) 

272 

273 def __setitem__(self, key, value): 

274 if not isinstance(value, VelocityProfile): 

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

276 self.profiles[key] = value 

277 

278 def __delitem__(self, key): 

279 self.profiles.remove(key) 

280 

281 def __getitem__(self, key): 

282 return self.profiles[key] 

283 

284 def __str__(self): 

285 rstr = "Container contains %d velocity profiles:\n\n" % self.nprofiles 

286 return rstr 

287 

288 @property 

289 def nprofiles(self): 

290 return len(self.profiles) 

291 

292 def append(self, value): 

293 if not isinstance(value, VelocityProfile): 

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

295 self.profiles.append(value) 

296 

297 def copy(self): 

298 return copy.deepcopy(self) 

299 

300 def lats(self): 

301 return num.array( 

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

303 

304 def lons(self): 

305 return num.array( 

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

307 

308 def _dataMatrix(self): 

309 if self.data_matrix is not None: 

310 return self.data_matrix 

311 

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

313 num.vstack([ 

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

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

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

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

318 ]), 

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

320 return self.data_matrix 

321 

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

323 ''' 

324 Create a regular sampled velocity matrix. 

325 

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

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

328 :type depth_range: tuple 

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

330 :type ddepth: float 

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

332 defaults to ``p`` 

333 :type phase: str 

334 :returns: Sample depths, veloctiy matrix 

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

336 ''' 

337 dmin, dmax = depth_range 

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

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

340 ndepth = sdepth.size 

341 

342 if uid not in self._velocity_matrix_cache: 

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

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

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

346 phase=phase) 

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

348 

349 return sdepth, self._velocity_matrix_cache[uid] 

350 

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

352 phase='p'): 

353 ''' 

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

355 

356 :param ref_profile: Reference profile 

357 :type ref_profile: :class:`VelocityProfile` 

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

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

360 :type depth_range: tuple, optional 

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

362 :type ddepth: float 

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

364 :type phase: str 

365 :returns: RMS factor length of N_profiles 

366 :rtype: :class:`numpy.ndarray` 

367 ''' 

368 if not isinstance(ref_profile, VelocityProfile): 

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

370 

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

372 phase=phase) 

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

374 

375 rms = num.empty(self.nprofiles) 

376 for p in range(self.nprofiles): 

377 profile = vel_matrix[p, :] 

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

379 return rms 

380 

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

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

383 ''' 

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

385 

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

387 

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

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

390 :type depth_range: tuple 

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

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

393 :type vel_range: tuple 

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

395 :type ddepth: float 

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

397 :type dvbin: float 

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

399 :type dvbin: float 

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

401 :type phase: str 

402 

403 :returns: :func:`numpy.histogram2d` 

404 :rtype: tuple 

405 ''' 

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

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

408 

409 # Velocity and depth bins 

410 if vel_range is None: 

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

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

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

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

415 

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

417 range=(vel_range, depth_range), 

418 bins=[nvbins, ndbins], 

419 normed=False) 

420 

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

422 ''' 

423 Get mean and standard deviation of velocity profile. 

424 

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

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

427 :type depth_range: tuple 

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

429 :type ddepth: float 

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

431 :type phase: str 

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

433 :rtype: tuple of :class:`numpy.ndarray` 

434 ''' 

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

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

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

438 

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

440 

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

442 ''' 

443 Get mode velocity profile and standard deviation. 

444 

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

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

447 :type depth_range: tuple 

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

449 :type ddepth: float 

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

451 :type phase: str 

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

453 :rtype: tuple of :class:`numpy.ndarray` 

454 ''' 

455 import scipy.stats 

456 

457 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth) 

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

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

460 

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

462 ''' 

463 Median velocity profile plus std variation. 

464 

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

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

467 :type depth_range: tuple 

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

469 :type ddepth: float 

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

471 :type phase: str 

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

473 :rtype: tuple of :class:`numpy.ndarray` 

474 ''' 

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

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

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

478 

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

480 

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

482 axes=None): 

483 ''' 

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

485 

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

487 :type vel_range: tuple, optional 

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

489 :type bins: int, optional 

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

491 defaults to 'vp' 

492 :type phase: str, optional 

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

494 :type figure: :class:`matplotlib.Figure`, optional 

495 ''' 

496 

497 import matplotlib.pyplot as plt 

498 

499 fig, ax = _getCanvas(axes) 

500 

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

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

503 

504 data = self._dataMatrix()[phase] 

505 

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

507 range=vel_range, bins=bins, 

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

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

510 transform=ax.transAxes, fontsize=10, 

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

512 

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

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

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

516 xscaled(1./km, ax) 

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

518 

519 if self.name is not None: 

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

521 

522 if axes is None: 

523 plt.show() 

524 

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

526 vel_range=None, dvbin=100., 

527 percent=False, 

528 plot_mode=True, plot_median=True, plot_mean=False, 

529 show_cbar=True, 

530 aspect=.02, 

531 phase='p', 

532 axes=None): 

533 ''' 

534 Plot a two 2D Histogram of seismic velocities. 

535 

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

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

538 :type depth_range: tuple 

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

540 :type vel_range: tuple 

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

542 :type ddepth: float 

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

544 :type dvbin: float 

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

546 :type dvbin: float 

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

548 :type phase: str 

549 :param plot_mode: Plot the Mode 

550 :type plot_mode: bool 

551 :param plot_mean: Plot the Mean 

552 :type plot_mean: bool 

553 :param plot_median: Plot the Median 

554 :type plot_median: bool 

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

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

557 ''' 

558 

559 import matplotlib.pyplot as plt 

560 

561 fig, ax = _getCanvas(axes) 

562 

563 ax = fig.gca() 

564 

565 if vel_range is not None: 

566 vmin, vmax = vel_range 

567 dmin, dmax = depth_range 

568 

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

570 depth_range=depth_range, 

571 ddepth=ddepth, dvbin=dvbin, 

572 ddbin=ddbin, phase=phase) 

573 vfield /= (ddbin / ddepth) 

574 

575 if percent: 

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

577 

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

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

580 interpolation='nearest', 

581 extent=grid_ext, aspect=aspect) 

582 

583 if show_cbar: 

584 cticks = num.unique( 

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

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

587 orientation='horizontal') 

588 if percent: 

589 cbar.set_label('Percent') 

590 else: 

591 cbar.set_label('Number of Profiles') 

592 

593 if plot_mode: 

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

595 ddepth=ddepth) 

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

597 sdepth[sdepth < dmax], 

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

599 

600 if plot_mean: 

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

602 ddepth=ddepth) 

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

604 sdepth[sdepth < dmax], 

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

606 

607 if plot_median: 

608 sdepth, vel_median, _ = self.medianVelocity( 

609 depth_range=depth_range, 

610 ddepth=ddepth) 

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

612 sdepth[sdepth < dmax], 

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

614 

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

616 

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

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

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

620 

621 ax.set_title('Crustal Velocity Distribution') 

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

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

624 yscaled(1./km, ax) 

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

626 ax.set_xlim(vel_range) 

627 

628 if self.name is not None: 

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

630 

631 if plot_mode or plot_mean or plot_median: 

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

633 leg.get_frame().set_alpha(.6) 

634 

635 if axes is None: 

636 plt.show() 

637 

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

639 ''' 

640 Plot a triangulated a depth surface exceeding velocity. 

641 ''' 

642 

643 import matplotlib.pyplot as plt 

644 

645 fig, ax = _getCanvas(axes) 

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

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

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

649 d = d[d > 0] 

650 

651 ax.tricontourf(lats, lons, d) 

652 

653 if axes is None: 

654 plt.show() 

655 

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

657 from pyrocko.plot import gmtpy 

658 lats = self.lats() 

659 lons = self.lons() 

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

661 

662 def darken(c, f=0.7): 

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

664 

665 gmt = gmtpy.GMT() 

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

667 J='M0/12', 

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

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

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

671 W='.2p') 

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

673 in_columns=[lons, lats], 

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

675 gmt.save(outfile) 

676 

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

678 ''' 

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

680 

681 :param v_max: maximal velocity 

682 :type vmax: float 

683 :param dz: depth is sampled in dz steps 

684 :type dz: float 

685 :param d_max: maximum depth 

686 :type d_max: int 

687 :param d_min: minimum depth 

688 :type d_min: int 

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

690 :rtype: list(num.array) 

691 ''' 

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

693 self.profile_exceed_velocity[:] = num.nan 

694 

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

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

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

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

699 continue 

700 if profile.vp[il] < v_max: 

701 continue 

702 else: 

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

704 break 

705 return self.profile_exceed_velocity 

706 

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

708 ''' 

709 Select profiles within a region by geographic corner coordinates. 

710 

711 :param west: west corner 

712 :type west: float 

713 :param east: east corner 

714 :type east: float 

715 :param south: south corner 

716 :type south: float 

717 :param north: north corner 

718 :type north: float 

719 :returns: Selected profiles 

720 :rtype: :class:`CrustDB` 

721 ''' 

722 r_container = self._emptyCopy() 

723 

724 for profile in self.profiles: 

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

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

727 r_container.append(profile) 

728 

729 return r_container 

730 

731 def selectPolygon(self, poly): 

732 ''' 

733 Select profiles within a polygon. 

734 

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

736 

737 :param poly: Latitude Longitude pairs of the polygon 

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

739 :returns: Selected profiles 

740 :rtype: :class:`CrustDB` 

741 ''' 

742 r_container = self._emptyCopy() 

743 

744 for profile in self.profiles: 

745 x = profile.lon 

746 y = profile.lat 

747 

748 inside = False 

749 p1x, p1y = poly[0] 

750 for p2x, p2y in poly: 

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

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

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

754 if p1y != p2y: 

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

756 (p2y - p1y) + p1x 

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

758 inside = not inside 

759 p1x, p1y = p2x, p2y 

760 if inside: 

761 r_container.append(profile) 

762 

763 return r_container 

764 

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

766 ''' 

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

768 

769 :param lat: Latitude in [deg] 

770 :type lat: float 

771 :param lon: Longitude in [deg] 

772 :type lon: float 

773 :param radius: Radius in [deg] 

774 :type radius: float 

775 :returns: Selected profiles 

776 :rtype: :class:`CrustDB` 

777 ''' 

778 r_container = self._emptyCopy() 

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

780 for profile in self.profiles: 

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

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

783 r_container.append(profile) 

784 

785 return r_container 

786 

787 def selectMinLayers(self, nlayers): 

788 ''' 

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

790 

791 :param nlayers: Minimum number of layers 

792 :type nlayers: int 

793 :returns: Selected profiles 

794 :rtype: :class:`CrustDB` 

795 ''' 

796 r_container = self._emptyCopy() 

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

798 

799 for profile in self.profiles: 

800 if profile.nlayers >= nlayers: 

801 r_container.append(profile) 

802 

803 return r_container 

804 

805 def selectMaxLayers(self, nlayers): 

806 ''' 

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

808 

809 :param nlayers: Maximum number of layers 

810 :type nlayers: int 

811 :returns: Selected profiles 

812 :rtype: :class:`CrustDB` 

813 ''' 

814 r_container = self._emptyCopy() 

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

816 

817 for profile in self.profiles: 

818 if profile.nlayers <= nlayers: 

819 r_container.append(profile) 

820 

821 return r_container 

822 

823 def selectMinDepth(self, depth): 

824 ''' 

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

826 

827 :param depth: Minumum depth in [m] 

828 :type depth: float 

829 :returns: Selected profiles 

830 :rtype: :class:`CrustDB` 

831 ''' 

832 r_container = self._emptyCopy() 

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

834 

835 for profile in self.profiles: 

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

837 r_container.append(profile) 

838 return r_container 

839 

840 def selectMaxDepth(self, depth): 

841 ''' 

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

843 

844 :param depth: Maximum depth in [m] 

845 :type depth: float 

846 :returns: Selected profiles 

847 :rtype: :class:`CrustDB` 

848 ''' 

849 r_container = self._emptyCopy() 

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

851 

852 for profile in self.profiles: 

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

854 r_container.append(profile) 

855 return r_container 

856 

857 def selectVp(self): 

858 ''' 

859 Select profiles describing P wave velocity. 

860 

861 :returns Selected profiles 

862 :rtype: :class:`CrustDB` 

863 ''' 

864 r_container = self._emptyCopy() 

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

866 

867 for profile in self.profiles: 

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

869 r_container.append(profile) 

870 return r_container 

871 

872 def selectVs(self): 

873 ''' 

874 Select profiles describing P wave velocity. 

875 

876 :returns: Selected profiles 

877 :rtype: :class:`CrustDB` 

878 ''' 

879 r_container = self._emptyCopy() 

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

881 

882 for profile in self.profiles: 

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

884 r_container.append(profile) 

885 return r_container 

886 

887 def _emptyCopy(self): 

888 r_container = CrustDB(parent=self) 

889 r_container.name = self.name 

890 return r_container 

891 

892 def exportCSV(self, filename=None): 

893 ''' 

894 Export as CSV file. 

895 

896 :param filename: Export filename 

897 :type filename: str 

898 ''' 

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

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

901 for profile in self.profiles: 

902 file.write(profile._csv()) 

903 

904 def exportYAML(self, filename=None): 

905 ''' 

906 Export as YAML file. 

907 

908 :param filename: Export filename 

909 :type filename: str 

910 ''' 

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

912 for profile in self.profiles: 

913 file.write(profile.__str__()) 

914 

915 @classmethod 

916 def readDatabase(cls, database_file): 

917 db = cls() 

918 CrustDB._read(db, database_file) 

919 return db 

920 

921 @staticmethod 

922 def _getRepositoryDatabase(): 

923 from pyrocko import config 

924 

925 name = path.basename(db_url) 

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

927 if not path.exists(data_path): 

928 from pyrocko import util 

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

930 

931 return data_path 

932 

933 def _read(self, database_file): 

934 ''' 

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

936 

937 File format: 

938 

939 uid lat/lon vp vs hc depth 

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

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

942 5.38 .00 12.50 7.30 c 

943 6.92 .00 13.20 19.80 c 

944 8.18 .00 .00 33.00 m 

945 

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

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

948 7.00 .00 18.50 19.50 

949 7.80 .00 .00 38.00 m 

950 

951 

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

953 

954 ''' 

955 

956 def get_empty_record(): 

957 meta = { 

958 'uid': num.nan, 

959 'geographical_location': None, 

960 'geological_province': None, 

961 'geological_age': None, 

962 'elevation': num.nan, 

963 'heatflow': num.nan, 

964 'measurement_method': None, 

965 'publication_reference': None 

966 } 

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

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

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

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

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

972 0., 0., meta) 

973 

974 nlayers = [] 

975 

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

977 if nlayer == 0: 

978 return 

979 self.append(VelocityProfile( 

980 vp=vp[:nlayer] * km, 

981 vs=vs[:nlayer] * km, 

982 h=h[:nlayer] * km, 

983 d=d[:nlayer] * km, 

984 lat=lat, lon=lon, 

985 **meta)) 

986 nlayers.append(nlayer) 

987 

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

989 ilayer = 0 

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

991 for line, dbline in enumerate(database): 

992 if dbline.isspace(): 

993 if not len(d) == 0: 

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

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

996 raise DatabaseError( 

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

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

999 

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

1001 ilayer = 0 

1002 else: 

1003 try: 

1004 if ilayer == 0: 

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

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

1007 lat = -lat 

1008 # Additional meta data 

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

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

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

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

1013 meta['heatflow'] = None 

1014 meta['geographical_location'] =\ 

1015 dbline[66:72].strip() 

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

1017 if ilayer == 1: 

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

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

1020 lon = -lon 

1021 # Additional meta data 

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

1023 meta['publication_reference'] =\ 

1024 dbline[66:72].strip() 

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

1026 try: 

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

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

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

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

1031 except ValueError: 

1032 pass 

1033 except ValueError: 

1034 logger.warning( 

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

1036 (line, dbline)) 

1037 while not database.readlines(): 

1038 pass 

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

1040 ilayer += 1 

1041 # Append last profile 

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

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

1044 (self.nprofiles, database_file))