1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5'''
6Access to the USGS Global Crustal Database.
8Simple queries and statistical analysis
9'''
11import numpy as num
12import copy
13import logging
14from os import path
16from pyrocko.guts import Object, String, Float, Int
17from pyrocko.guts_array import Array
19from pyrocko.cake import LayeredModel, Material
20from pyrocko.plot.cake_plot import my_model_plot, xscaled, yscaled
22from .crustdb_abbr import ageKey, provinceKey, referenceKey, pubYear # noqa
24logger = logging.getLogger('pyrocko.dataset.crustdb')
25THICKNESS_HALFSPACE = 2
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}
37class DatabaseError(Exception):
38 pass
41class ProfileEmpty(Exception):
42 pass
45def _getCanvas(axes):
47 import matplotlib.pyplot as plt
49 if axes is None:
50 fig = plt.figure()
51 return fig, fig.gca()
52 return axes.figure, axes
55def xoffset_scale(offset, scale, ax):
56 from matplotlib.ticker import ScalarFormatter, AutoLocator
58 class FormatVelocities(ScalarFormatter):
59 @staticmethod
60 def __call__(value, pos):
61 return u'%.1f' % ((value-offset) * scale)
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)]
68 ax.get_xaxis().set_major_formatter(FormatVelocities())
69 ax.get_xaxis().set_major_locator(OffsetLocator())
72class VelocityProfile(Object):
73 uid = Int.T(
74 optional=True,
75 help='Unique ID of measurement')
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]')
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')
118 def __init__(self, *args, **kwargs):
119 Object.__init__(self, *args, **kwargs)
121 self.h = num.abs(self.d - num.roll(self.d, -1))
122 self.h[-1] = 0
123 self.nlayers = self.h.size
125 self.geographical_location = '%s (%s)' % (
126 provinceKey(self.geographical_location),
127 self.geographical_location)
129 self.vs[self.vs == 0] = num.nan
130 self.vp[self.vp == 0] = num.nan
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
137 @property
138 def publication_year__(self):
139 return pubYear(self.publication_reference)
141 def interpolateProfile(self, depths, phase='p', stepped=True):
142 '''
143 Get a continuous velocity function at arbitrary depth.
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 '''
155 if phase not in ['s', 'p']:
156 raise AttributeError('Phase has to be either \'p\' or \'s\'.')
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
164 if vel.size == 0:
165 raise ProfileEmpty('Phase %s does not contain velocities' % phase)
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.')
173 return res
175 def plot(self, axes=None):
176 '''
177 Plot the velocity profile, see :class:`pyrocko.cake`.
179 :param axes: Axes to plot into.
180 :type axes: :class:`matplotlib.Axes`
181 '''
183 import matplotlib.pyplot as plt
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()
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, ''
201 return LayeredModel.from_scanlines(iterLines())
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])
211 @property
212 def geog_loc_long(self):
213 return provinceKey(self.geog_loc)
215 @property
216 def geol_age_long(self):
217 return ageKey(self.geol_age)
219 @property
220 def has_s(self):
221 return num.any(self.vp)
223 @property
224 def has_p(self):
225 return num.any(self.vs)
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
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
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 '''
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
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())
270 def __len__(self):
271 return len(self.profiles)
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
278 def __delitem__(self, key):
279 self.profiles.remove(key)
281 def __getitem__(self, key):
282 return self.profiles[key]
284 def __str__(self):
285 rstr = "Container contains %d velocity profiles:\n\n" % self.nprofiles
286 return rstr
288 @property
289 def nprofiles(self):
290 return len(self.profiles)
292 def append(self, value):
293 if not isinstance(value, VelocityProfile):
294 raise TypeError('Element is not a VelocityProfile')
295 self.profiles.append(value)
297 def copy(self):
298 return copy.deepcopy(self)
300 def lats(self):
301 return num.array(
302 [p.lat for p in self.profiles])
304 def lons(self):
305 return num.array(
306 [p.lon for p in self.profiles])
308 def _dataMatrix(self):
309 if self.data_matrix is not None:
310 return self.data_matrix
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
322 def velocityMatrix(self, depth_range=(0, 60000.), ddepth=100., phase='p'):
323 '''
324 Create a regular sampled velocity matrix.
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
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)
349 return sdepth, self._velocity_matrix_cache[uid]
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.
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')
371 sdepth, vel_matrix = self.velocityMatrix(depth_range, ddepth,
372 phase=phase)
373 ref_vel = ref_profile.interpolateProfile(sdepth, phase=phase)
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
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.
386 Check :func:`numpy.histogram2d` for more information.
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
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)
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)
416 return num.histogram2d(v_mat.flatten(), d_vec,
417 range=(vel_range, depth_range),
418 bins=[nvbins, ndbins],
419 normed=False)
421 def meanVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'):
422 '''
423 Get mean and standard deviation of velocity profile.
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)
439 return sdepth, v_mean.flatten(), v_std.flatten()
441 def modeVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'):
442 '''
443 Get mode velocity profile and standard deviation.
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
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()
461 def medianVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'):
462 '''
463 Median velocity profile plus std variation.
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)
479 return sdepth, v_mean.flatten(), v_std.flatten()
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.
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 '''
497 import matplotlib.pyplot as plt
499 fig, ax = _getCanvas(axes)
501 if phase not in ['vp', 'vs']:
502 raise AttributeError('phase has to be either vp or vs')
504 data = self._dataMatrix()[phase]
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)
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)
519 if self.name is not None:
520 ax.set_title('%s for %s' % (ax.get_title(), self.name))
522 if axes is None:
523 plt.show()
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.
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 '''
559 import matplotlib.pyplot as plt
561 fig, ax = _getCanvas(axes)
563 ax = fig.gca()
565 if vel_range is not None:
566 vmin, vmax = vel_range
567 dmin, dmax = depth_range
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)
575 if percent:
576 vfield /= vfield.sum(axis=1)[num.newaxis, :]
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)
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')
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')
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')
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')
615 ax.grid(True, which="both", color="w", linewidth=.8, alpha=.4)
617 ax.text(.025, .025, '%d Profiles' % self.nprofiles,
618 color='w', alpha=.7,
619 transform=ax.transAxes, fontsize=9, va='bottom', ha='left')
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)
628 if self.name is not None:
629 ax.set_title('%s for %s' % (ax.get_title(), self.name))
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)
635 if axes is None:
636 plt.show()
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 '''
643 import matplotlib.pyplot as plt
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]
651 ax.tricontourf(lats, lons, d)
653 if axes is None:
654 plt.show()
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())
662 def darken(c, f=0.7):
663 return (c[0]*f, c[1]*f, c[2]*f)
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)
677 def exceedVelocity(self, v_max, d_min=0, d_max=60):
678 '''
679 Returns the last depth ``v_max`` has not been exceeded.
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
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
707 def selectRegion(self, west, east, south, north):
708 '''
709 Select profiles within a region by geographic corner coordinates.
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()
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)
729 return r_container
731 def selectPolygon(self, poly):
732 '''
733 Select profiles within a polygon.
735 The algorithm is called the **Ray Casting Method**
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()
744 for profile in self.profiles:
745 x = profile.lon
746 y = profile.lat
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)
763 return r_container
765 def selectLocation(self, lat, lon, radius=10):
766 '''
767 Select profiles at a geographic location within a ``radius``.
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)
785 return r_container
787 def selectMinLayers(self, nlayers):
788 '''
789 Select profiles with more than ``nlayers``.
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)
799 for profile in self.profiles:
800 if profile.nlayers >= nlayers:
801 r_container.append(profile)
803 return r_container
805 def selectMaxLayers(self, nlayers):
806 '''
807 Select profiles with more than ``nlayers``.
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)
817 for profile in self.profiles:
818 if profile.nlayers <= nlayers:
819 r_container.append(profile)
821 return r_container
823 def selectMinDepth(self, depth):
824 '''
825 Select profiles describing layers deeper than ``depth``.
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)
835 for profile in self.profiles:
836 if profile.d.max() >= depth:
837 r_container.append(profile)
838 return r_container
840 def selectMaxDepth(self, depth):
841 '''
842 Select profiles describing layers shallower than ``depth``.
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)
852 for profile in self.profiles:
853 if profile.d.max() <= depth:
854 r_container.append(profile)
855 return r_container
857 def selectVp(self):
858 '''
859 Select profiles describing P wave velocity.
861 :returns Selected profiles
862 :rtype: :class:`CrustDB`
863 '''
864 r_container = self._emptyCopy()
865 logger.info('Selecting profiles providing Vp...')
867 for profile in self.profiles:
868 if not num.all(num.isnan(profile.vp)):
869 r_container.append(profile)
870 return r_container
872 def selectVs(self):
873 '''
874 Select profiles describing P wave velocity.
876 :returns: Selected profiles
877 :rtype: :class:`CrustDB`
878 '''
879 r_container = self._emptyCopy()
880 logger.info('Selecting profiles providing Vs...')
882 for profile in self.profiles:
883 if not num.all(num.isnan(profile.vs)):
884 r_container.append(profile)
885 return r_container
887 def _emptyCopy(self):
888 r_container = CrustDB(parent=self)
889 r_container.name = self.name
890 return r_container
892 def exportCSV(self, filename=None):
893 '''
894 Export as CSV file.
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())
904 def exportYAML(self, filename=None):
905 '''
906 Export as YAML file.
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__())
915 @classmethod
916 def readDatabase(cls, database_file):
917 db = cls()
918 CrustDB._read(db, database_file)
919 return db
921 @staticmethod
922 def _getRepositoryDatabase():
923 from pyrocko import config
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)
931 return data_path
933 def _read(self, database_file):
934 '''
935 Reads in the the GSN database file and puts it in CrustDB.
937 File format:
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
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
952 :param database_file: path to database file, type string
954 '''
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)
974 nlayers = []
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)
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)
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))