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'''
10from __future__ import absolute_import
12import numpy as num
13import copy
14import logging
15from os import path
17from pyrocko.guts import Object, String, Float, Int
18from pyrocko.guts_array import Array
20from pyrocko.cake import LayeredModel, Material
21from pyrocko.plot.cake_plot import my_model_plot, xscaled, yscaled
23from .crustdb_abbr import ageKey, provinceKey, referenceKey, pubYear # noqa
25logger = logging.getLogger('pyrocko.dataset.crustdb')
26THICKNESS_HALFSPACE = 2
28db_url = 'https://mirror.pyrocko.org/gsc20130501.txt'
29km = 1e3
30vel_labels = {
31 'vp': '$V_P$',
32 'p': '$V_P$',
33 'vs': '$V_S$',
34 's': '$V_S$',
35}
38class DatabaseError(Exception):
39 pass
42class ProfileEmpty(Exception):
43 pass
46def _getCanvas(axes):
48 import matplotlib.pyplot as plt
50 if axes is None:
51 fig = plt.figure()
52 return fig, fig.gca()
53 return axes.figure, axes
56def xoffset_scale(offset, scale, ax):
57 from matplotlib.ticker import ScalarFormatter, AutoLocator
59 class FormatVelocities(ScalarFormatter):
60 @staticmethod
61 def __call__(value, pos):
62 return u'%.1f' % ((value-offset) * scale)
64 class OffsetLocator(AutoLocator):
65 def tick_values(self, vmin, vmax):
66 return [v + offset for v in
67 AutoLocator.tick_values(self, vmin, vmax)]
69 ax.get_xaxis().set_major_formatter(FormatVelocities())
70 ax.get_xaxis().set_major_locator(OffsetLocator())
73class VelocityProfile(Object):
74 uid = Int.T(
75 optional=True,
76 help='Unique ID of measurement')
78 lat = Float.T(
79 help='Latitude [deg]')
80 lon = Float.T(
81 help='Longitude [deg]')
82 elevation = Float.T(
83 default=num.nan,
84 help='Elevation [m]')
85 vp = Array.T(
86 shape=(None, 1),
87 help='P Wave velocities [m/s]')
88 vs = Array.T(
89 shape=(None, 1),
90 help='S Wave velocities [m/s]')
91 d = Array.T(
92 shape=(None, 1),
93 help='Interface depth, top [m]')
94 h = Array.T(
95 shape=(None, 1),
96 help='Interface thickness [m]')
98 heatflow = Float.T(
99 optional=True,
100 help='Heatflow [W/m^2]')
101 geographical_location = String.T(
102 optional=True,
103 help='Geographic Location')
104 geological_province = String.T(
105 optional=True,
106 help='Geological Province')
107 geological_age = String.T(
108 optional=True,
109 help='Geological Age')
110 measurement_method = Int.T(
111 optional=True,
112 help='Measurement method')
113 publication_reference = String.T(
114 optional=True,
115 help='Publication Reference')
116 publication_year__ = Int.T(
117 help='Publication Date')
119 def __init__(self, *args, **kwargs):
120 Object.__init__(self, *args, **kwargs)
122 self.h = num.abs(self.d - num.roll(self.d, -1))
123 self.h[-1] = 0
124 self.nlayers = self.h.size
126 self.geographical_location = '%s (%s)' % (
127 provinceKey(self.geographical_location),
128 self.geographical_location)
130 self.vs[self.vs == 0] = num.nan
131 self.vp[self.vp == 0] = num.nan
133 self._step_vp = num.repeat(self.vp, 2)
134 self._step_vs = num.repeat(self.vs, 2)
135 self._step_d = num.roll(num.repeat(self.d, 2), -1)
136 self._step_d[-1] = self._step_d[-2] + THICKNESS_HALFSPACE
138 @property
139 def publication_year__(self):
140 return pubYear(self.publication_reference)
142 def interpolateProfile(self, depths, phase='p', stepped=True):
143 '''
144 Get a continuous velocity function at arbitrary depth.
146 :param depth: Depths to interpolate
147 :type depth: :class:`numpy.ndarray`
148 :param phase: P or S wave velocity, **p** or **s**
149 :type phase: str, optional
150 :param stepped: Use a stepped velocity function or gradient
151 :type stepped: bool
152 :returns: velocities at requested depths
153 :rtype: :class:`numpy.ndarray`
154 '''
156 if phase not in ['s', 'p']:
157 raise AttributeError('Phase has to be either \'p\' or \'s\'.')
159 if phase == 'p':
160 vel = self._step_vp if stepped else self.vp
161 elif phase == 's':
162 vel = self._step_vs if stepped else self.vs
163 d = self._step_d if stepped else self.d
165 if vel.size == 0:
166 raise ProfileEmpty('Phase %s does not contain velocities' % phase)
168 try:
169 res = num.interp(depths, d, vel,
170 left=num.nan, right=num.nan)
171 except ValueError:
172 raise ValueError('Could not interpolate velocity profile.')
174 return res
176 def plot(self, axes=None):
177 '''
178 Plot the velocity profile, see :class:`pyrocko.cake`.
180 :param axes: Axes to plot into.
181 :type axes: :class:`matplotlib.Axes`
182 '''
184 import matplotlib.pyplot as plt
186 fig, ax = _getCanvas(axes)
187 my_model_plot(self.getLayeredModel(), axes=axes)
188 ax.set_title('Global Crustal Database\n'
189 'Velocity Structure at {p.lat:.4f}N, '
190 ' {p.lat:.4f}E (uid {p.uid})'.format(p=self))
191 if axes is None:
192 plt.show()
194 def getLayeredModel(self):
195 '''
196 Get a layered model, see :class:`pyrocko.cake.LayeredModel`.
197 '''
198 def iterLines():
199 for il, m in enumerate(self.iterLayers()):
200 yield self.d[il], m, ''
202 return LayeredModel.from_scanlines(iterLines())
204 def iterLayers(self):
205 '''
206 Iterator yielding a :class:`pyrocko.cake.Material` for each layer.
207 '''
208 for il in range(self.nlayers):
209 yield Material(vp=self.vp[il],
210 vs=self.vs[il])
212 @property
213 def geog_loc_long(self):
214 return provinceKey(self.geog_loc)
216 @property
217 def geol_age_long(self):
218 return ageKey(self.geol_age)
220 @property
221 def has_s(self):
222 return num.any(self.vp)
224 @property
225 def has_p(self):
226 return num.any(self.vs)
228 def get_weeded(self):
229 '''
230 Get weeded representation of layers used in the profile.
231 See :func:`pyrocko.cake.get_weeded` for details.
232 '''
233 weeded = num.zeros((self.nlayers, 4))
234 weeded[:, 0] = self.d
235 weeded[:, 1] = self.vp
236 weeded[:, 2] = self.vs
238 def _csv(self):
239 output = ''
240 for d in range(len(self.h)):
241 output += (
242 '{p.uid}, {p.lat}, {p.lon},'
243 ' {vp}, {vs}, {h}, {d}, {p.publication_reference}\n'
244 ).format(
245 p=self,
246 vp=self.vp[d], vs=self.vs[d], h=self.h[d], d=self.d[d])
247 return output
250class CrustDB(object):
251 '''
252 CrustDB is a container for :class:`VelocityProfile` and provides
253 functions for spatial selection, querying, processing and visualising
254 data from the Global Crustal Database.
255 '''
257 def __init__(self, database_file=None, parent=None):
258 self.profiles = []
259 self._velocity_matrix_cache = {}
260 self.data_matrix = None
261 self.name = None
262 self.database_file = database_file
264 if parent is not None:
265 pass
266 elif database_file is not None:
267 self._read(database_file)
268 else:
269 self._read(self._getRepositoryDatabase())
271 def __len__(self):
272 return len(self.profiles)
274 def __setitem__(self, key, value):
275 if not isinstance(value, VelocityProfile):
276 raise TypeError('Element is not a VelocityProfile')
277 self.profiles[key] = value
279 def __delitem__(self, key):
280 self.profiles.remove(key)
282 def __getitem__(self, key):
283 return self.profiles[key]
285 def __str__(self):
286 rstr = "Container contains %d velocity profiles:\n\n" % self.nprofiles
287 return rstr
289 @property
290 def nprofiles(self):
291 return len(self.profiles)
293 def append(self, value):
294 if not isinstance(value, VelocityProfile):
295 raise TypeError('Element is not a VelocityProfile')
296 self.profiles.append(value)
298 def copy(self):
299 return copy.deepcopy(self)
301 def lats(self):
302 return num.array(
303 [p.lat for p in self.profiles])
305 def lons(self):
306 return num.array(
307 [p.lon for p in self.profiles])
309 def _dataMatrix(self):
310 if self.data_matrix is not None:
311 return self.data_matrix
313 self.data_matrix = num.core.records.fromarrays(
314 num.vstack([
315 num.concatenate([p.vp for p in self.profiles]),
316 num.concatenate([p.vs for p in self.profiles]),
317 num.concatenate([p.h for p in self.profiles]),
318 num.concatenate([p.d for p in self.profiles])
319 ]),
320 names='vp, vs, h, d')
321 return self.data_matrix
323 def velocityMatrix(self, depth_range=(0, 60000.), ddepth=100., phase='p'):
324 '''
325 Create a regular sampled velocity matrix.
327 :param depth_range: Depth range, ``(dmin, dmax)``,
328 defaults to ``(0, 6000.)``
329 :type depth_range: tuple
330 :param ddepth: Stepping in [m], defaults to ``100.``
331 :type ddepth: float
332 :param phase: Phase to calculate ``p`` or ``s``,
333 defaults to ``p``
334 :type phase: str
335 :returns: Sample depths, veloctiy matrix
336 :rtype: tuple, (sample_depth, :class:`numpy.ndarray`)
337 '''
338 dmin, dmax = depth_range
339 uid = '.'.join(map(repr, (dmin, dmax, ddepth, phase)))
340 sdepth = num.linspace(dmin, dmax, (dmax - dmin) / ddepth)
341 ndepth = sdepth.size
343 if uid not in self._velocity_matrix_cache:
344 vel_mat = num.empty((self.nprofiles, ndepth))
345 for ip, profile in enumerate(self.profiles):
346 vel_mat[ip, :] = profile.interpolateProfile(sdepth,
347 phase=phase)
348 self._velocity_matrix_cache[uid] = num.ma.masked_invalid(vel_mat)
350 return sdepth, self._velocity_matrix_cache[uid]
352 def rmsRank(self, ref_profile, depth_range=(0, 3500.), ddepth=100.,
353 phase='p'):
354 '''
355 Correlates ``ref_profile`` to each profile in the database.
357 :param ref_profile: Reference profile
358 :type ref_profile: :class:`VelocityProfile`
359 :param depth_range: Depth range in [m], ``(dmin, dmax)``,
360 defaults to ``(0, 35000.)``
361 :type depth_range: tuple, optional
362 :param ddepth: Stepping in [m], defaults to ``100.``
363 :type ddepth: float
364 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p``
365 :type phase: str
366 :returns: RMS factor length of N_profiles
367 :rtype: :class:`numpy.ndarray`
368 '''
369 if not isinstance(ref_profile, VelocityProfile):
370 raise ValueError('ref_profile is not a VelocityProfile')
372 sdepth, vel_matrix = self.velocityMatrix(depth_range, ddepth,
373 phase=phase)
374 ref_vel = ref_profile.interpolateProfile(sdepth, phase=phase)
376 rms = num.empty(self.nprofiles)
377 for p in range(self.nprofiles):
378 profile = vel_matrix[p, :]
379 rms[p] = num.sqrt(profile**2 - ref_vel**2).sum() / ref_vel.size
380 return rms
382 def histogram2d(self, depth_range=(0., 60000.), vel_range=None,
383 ddepth=100., dvbin=100., ddbin=2000., phase='p'):
384 '''
385 Create a 2D Histogram of all the velocity profiles.
387 Check :func:`numpy.histogram2d` for more information.
389 :param depth_range: Depth range in [m], ``(dmin, dmax)``,
390 defaults to ``(0., 60000.)``
391 :type depth_range: tuple
392 :param vel_range: Depth range, ``(vmin, vmax)``,
393 defaults to ``(5500., 8500.)``
394 :type vel_range: tuple
395 :param ddepth: Stepping in [km], defaults to ``100.``
396 :type ddepth: float
397 :param dvbin: Bin size in velocity dimension [m/s], defaults to 100.
398 :type dvbin: float
399 :param dvbin: Bin size in depth dimension [m], defaults to 2000.
400 :type dvbin: float
401 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p``
402 :type phase: str
404 :returns: :func:`numpy.histogram2d`
405 :rtype: tuple
406 '''
407 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth, phase=phase)
408 d_vec = num.tile(sdepth, self.nprofiles)
410 # Velocity and depth bins
411 if vel_range is None:
412 vel_range = ((v_mat.min() // 1e2) * 1e2,
413 (v_mat.max() // 1e2) * 1e2)
414 nvbins = int((vel_range[1] - vel_range[0]) / dvbin)
415 ndbins = int((depth_range[1] - depth_range[0]) / ddbin)
417 return num.histogram2d(v_mat.flatten(), d_vec,
418 range=(vel_range, depth_range),
419 bins=[nvbins, ndbins],
420 normed=False)
422 def meanVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'):
423 '''
424 Get mean and standard deviation of velocity profile.
426 :param depth_range: Depth range in [m], ``(dmin, dmax)``,
427 defaults to ``(0., 60000.)``
428 :type depth_range: tuple
429 :param ddepth: Stepping in [m], defaults to ``100.``
430 :type ddepth: float
431 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p``
432 :type phase: str
433 :returns: depth vector, mean velocities, standard deviations
434 :rtype: tuple of :class:`numpy.ndarray`
435 '''
436 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth, phase=phase)
437 v_mean = num.ma.mean(v_mat, axis=0)
438 v_std = num.ma.std(v_mat, axis=0)
440 return sdepth, v_mean.flatten(), v_std.flatten()
442 def modeVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'):
443 '''
444 Get mode velocity profile and standard deviation.
446 :param depth_range: Depth range in [m], ``(dmin, dmax)``,
447 defaults to ``(0., 60000.)``
448 :type depth_range: tuple
449 :param ddepth: Stepping in [m], defaults to ``100.``
450 :type ddepth: float
451 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p``
452 :type phase: str
453 :returns: depth vector, mode velocity, number of counts at each depth
454 :rtype: tuple of :class:`numpy.ndarray`
455 '''
456 import scipy.stats
458 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth)
459 v_mode, v_counts = scipy.stats.mstats.mode(v_mat, axis=0)
460 return sdepth, v_mode.flatten(), v_counts.flatten()
462 def medianVelocity(self, depth_range=(0., 60000.), ddepth=100., phase='p'):
463 '''
464 Median velocity profile plus std variation.
466 :param depth_range: Depth range in [m], ``(dmin, dmax)``,
467 defaults to ``(0., 60000.)``
468 :type depth_range: tuple
469 :param ddepth: Stepping in [m], defaults to ``100.``
470 :type ddepth: float
471 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p``
472 :type phase: str
473 :returns: depth vector, median velocities, standard deviations
474 :rtype: tuple of :class:`numpy.ndarray`
475 '''
476 sdepth, v_mat = self.velocityMatrix(depth_range, ddepth, phase=phase)
477 v_mean = num.ma.median(v_mat, axis=0)
478 v_std = num.ma.std(v_mat, axis=0)
480 return sdepth, v_mean.flatten(), v_std.flatten()
482 def plotHistogram(self, vel_range=None, bins=36, phase='vp',
483 axes=None):
484 '''
485 Plot 1D histogram of seismic velocities in the container.
487 :param vel_range: Velocity range, defaults to (5.5, 8.5)
488 :type vel_range: tuple, optional
489 :param bins: bins, defaults to 30 (see :func:`numpy.histogram`)
490 :type bins: int, optional
491 :param phase: Property to plot out of ``['vp', 'vs']``,
492 defaults to 'vp'
493 :type phase: str, optional
494 :param figure: Figure to plot in, defaults to None
495 :type figure: :class:`matplotlib.Figure`, optional
496 '''
498 import matplotlib.pyplot as plt
500 fig, ax = _getCanvas(axes)
502 if phase not in ['vp', 'vs']:
503 raise AttributeError('phase has to be either vp or vs')
505 data = self._dataMatrix()[phase]
507 ax.hist(data, weights=self.data_matrix['h'],
508 range=vel_range, bins=bins,
509 color='g', alpha=.5)
510 ax.text(.95, .95, '%d Profiles' % self.nprofiles,
511 transform=ax.transAxes, fontsize=10,
512 va='top', ha='right', alpha=.7)
514 ax.set_title('Distribution of %s' % vel_labels[phase])
515 ax.set_xlabel('%s [km/s]' % vel_labels[phase])
516 ax.set_ylabel('Cumulative occurrence [N]')
517 xscaled(1./km, ax)
518 ax.yaxis.grid(alpha=.4)
520 if self.name is not None:
521 ax.set_title('%s for %s' % (ax.get_title(), self.name))
523 if axes is None:
524 plt.show()
526 def plot(self, depth_range=(0, 60000.), ddepth=100., ddbin=2000.,
527 vel_range=None, dvbin=100.,
528 percent=False,
529 plot_mode=True, plot_median=True, plot_mean=False,
530 show_cbar=True,
531 aspect=.02,
532 phase='p',
533 axes=None):
534 '''
535 Plot a two 2D Histogram of seismic velocities.
537 :param depth_range: Depth range, ``(dmin, dmax)``,
538 defaults to ``(0, 60)``
539 :type depth_range: tuple
540 :param vel_range: Velocity range, ``(vmin, vmax)``
541 :type vel_range: tuple
542 :param ddepth: Stepping in [m], defaults to ``.1``
543 :type ddepth: float
544 :param dvbin: Bin size in velocity dimension [m/s], defaults to .1
545 :type dvbin: float
546 :param dvbin: Bin size in depth dimension [m], defaults to 2000.
547 :type dvbin: float
548 :param phase: Phase to calculate ``p`` or ``s``, defaults to ``p``
549 :type phase: str
550 :param plot_mode: Plot the Mode
551 :type plot_mode: bool
552 :param plot_mean: Plot the Mean
553 :type plot_mean: bool
554 :param plot_median: Plot the Median
555 :type plot_median: bool
556 :param axes: Axes to plot into, defaults to None
557 :type axes: :class:`matplotlib.Axes`
558 '''
560 import matplotlib.pyplot as plt
562 fig, ax = _getCanvas(axes)
564 ax = fig.gca()
566 if vel_range is not None:
567 vmin, vmax = vel_range
568 dmin, dmax = depth_range
570 vfield, vedg, dedg = self.histogram2d(vel_range=vel_range,
571 depth_range=depth_range,
572 ddepth=ddepth, dvbin=dvbin,
573 ddbin=ddbin, phase=phase)
574 vfield /= (ddbin / ddepth)
576 if percent:
577 vfield /= vfield.sum(axis=1)[num.newaxis, :]
579 grid_ext = [vedg[0], vedg[-1], dedg[-1], dedg[0]]
580 histogram = ax.imshow(vfield.swapaxes(0, 1),
581 interpolation='nearest',
582 extent=grid_ext, aspect=aspect)
584 if show_cbar:
585 cticks = num.unique(
586 num.arange(0, vfield.max(), vfield.max() // 10).round())
587 cbar = fig.colorbar(histogram, ticks=cticks, format='%1i',
588 orientation='horizontal')
589 if percent:
590 cbar.set_label('Percent')
591 else:
592 cbar.set_label('Number of Profiles')
594 if plot_mode:
595 sdepth, vel_mode, _ = self.modeVelocity(depth_range=depth_range,
596 ddepth=ddepth)
597 ax.plot(vel_mode[sdepth < dmax] + ddepth/2,
598 sdepth[sdepth < dmax],
599 alpha=.8, color='w', label='Mode')
601 if plot_mean:
602 sdepth, vel_mean, _ = self.meanVelocity(depth_range=depth_range,
603 ddepth=ddepth)
604 ax.plot(vel_mean[sdepth < dmax] + ddepth/2,
605 sdepth[sdepth < dmax],
606 alpha=.8, color='w', linestyle='--', label='Mean')
608 if plot_median:
609 sdepth, vel_median, _ = self.medianVelocity(
610 depth_range=depth_range,
611 ddepth=ddepth)
612 ax.plot(vel_median[sdepth < dmax] + ddepth/2,
613 sdepth[sdepth < dmax],
614 alpha=.8, color='w', linestyle=':', label='Median')
616 ax.grid(True, which="both", color="w", linewidth=.8, alpha=.4)
618 ax.text(.025, .025, '%d Profiles' % self.nprofiles,
619 color='w', alpha=.7,
620 transform=ax.transAxes, fontsize=9, va='bottom', ha='left')
622 ax.set_title('Crustal Velocity Distribution')
623 ax.set_xlabel('%s [km/s]' % vel_labels[phase])
624 ax.set_ylabel('Depth [km]')
625 yscaled(1./km, ax)
626 xoffset_scale(dvbin/2, 1./km, ax)
627 ax.set_xlim(vel_range)
629 if self.name is not None:
630 ax.set_title('%s for %s' % (ax.get_title(), self.name))
632 if plot_mode or plot_mean or plot_median:
633 leg = ax.legend(loc=1, fancybox=True, prop={'size': 10.})
634 leg.get_frame().set_alpha(.6)
636 if axes is None:
637 plt.show()
639 def plotVelocitySurface(self, v_max, d_min=0., d_max=6000., axes=None):
640 '''
641 Plot a triangulated a depth surface exceeding velocity.
642 '''
644 import matplotlib.pyplot as plt
646 fig, ax = _getCanvas(axes)
647 d = self.exceedVelocity(v_max, d_min, d_max)
648 lons = self.lons()[d > 0]
649 lats = self.lats()[d > 0]
650 d = d[d > 0]
652 ax.tricontourf(lats, lons, d)
654 if axes is None:
655 plt.show()
657 def plotMap(self, outfile, **kwargs):
658 from pyrocko.plot import gmtpy
659 lats = self.lats()
660 lons = self.lons()
661 s, n, w, e = (lats.min(), lats.max(), lons.min(), lons.max())
663 def darken(c, f=0.7):
664 return (c[0]*f, c[1]*f, c[2]*f)
666 gmt = gmtpy.GMT()
667 gmt.psbasemap(B='40/20',
668 J='M0/12',
669 R='%f/%f/%f/%f' % (w, e, s, n))
670 gmt.pscoast(R=True, J=True,
671 D='i', S='216/242/254', A=10000,
672 W='.2p')
673 gmt.psxy(R=True, J=True,
674 in_columns=[lons, lats],
675 S='c2p', G='black')
676 gmt.save(outfile)
678 def exceedVelocity(self, v_max, d_min=0, d_max=60):
679 '''
680 Returns the last depth ``v_max`` has not been exceeded.
682 :param v_max: maximal velocity
683 :type vmax: float
684 :param dz: depth is sampled in dz steps
685 :type dz: float
686 :param d_max: maximum depth
687 :type d_max: int
688 :param d_min: minimum depth
689 :type d_min: int
690 :returns: Lat, Lon, Depth and uid where ``v_max`` is exceeded
691 :rtype: list(num.array)
692 '''
693 self.profile_exceed_velocity = num.empty(len(self.profiles))
694 self.profile_exceed_velocity[:] = num.nan
696 for ip, profile in enumerate(self.profiles):
697 for il in range(len(profile.d)):
698 if profile.d[il] <= d_min\
699 or profile.d[il] >= d_max:
700 continue
701 if profile.vp[il] < v_max:
702 continue
703 else:
704 self.profile_exceed_velocity[ip] = profile.d[il]
705 break
706 return self.profile_exceed_velocity
708 def selectRegion(self, west, east, south, north):
709 '''
710 Select profiles within a region by geographic corner coordinates.
712 :param west: west corner
713 :type west: float
714 :param east: east corner
715 :type east: float
716 :param south: south corner
717 :type south: float
718 :param north: north corner
719 :type north: float
720 :returns: Selected profiles
721 :rtype: :class:`CrustDB`
722 '''
723 r_container = self._emptyCopy()
725 for profile in self.profiles:
726 if profile.lon >= west and profile.lon <= east \
727 and profile.lat <= north and profile.lat >= south:
728 r_container.append(profile)
730 return r_container
732 def selectPolygon(self, poly):
733 '''
734 Select profiles within a polygon.
736 The algorithm is called the **Ray Casting Method**
738 :param poly: Latitude Longitude pairs of the polygon
739 :type param: list of :class:`numpy.ndarray`
740 :returns: Selected profiles
741 :rtype: :class:`CrustDB`
742 '''
743 r_container = self._emptyCopy()
745 for profile in self.profiles:
746 x = profile.lon
747 y = profile.lat
749 inside = False
750 p1x, p1y = poly[0]
751 for p2x, p2y in poly:
752 if y >= min(p1y, p2y):
753 if y <= max(p1y, p2y):
754 if x <= max(p1x, p2x):
755 if p1y != p2y:
756 xints = (y - p1y) * (p2x - p1x) / \
757 (p2y - p1y) + p1x
758 if p1x == p2x or x <= xints:
759 inside = not inside
760 p1x, p1y = p2x, p2y
761 if inside:
762 r_container.append(profile)
764 return r_container
766 def selectLocation(self, lat, lon, radius=10):
767 '''
768 Select profiles at a geographic location within a ``radius``.
770 :param lat: Latitude in [deg]
771 :type lat: float
772 :param lon: Longitude in [deg]
773 :type lon: float
774 :param radius: Radius in [deg]
775 :type radius: float
776 :returns: Selected profiles
777 :rtype: :class:`CrustDB`
778 '''
779 r_container = self._emptyCopy()
780 logger.info('Selecting location %f, %f (r=%f)...' % (lat, lon, radius))
781 for profile in self.profiles:
782 if num.sqrt((lat - profile.lat)**2 +
783 (lon - profile.lon)**2) <= radius:
784 r_container.append(profile)
786 return r_container
788 def selectMinLayers(self, nlayers):
789 '''
790 Select profiles with more than ``nlayers``.
792 :param nlayers: Minimum number of layers
793 :type nlayers: int
794 :returns: Selected profiles
795 :rtype: :class:`CrustDB`
796 '''
797 r_container = self._emptyCopy()
798 logger.info('Selecting minimum %d layers...' % nlayers)
800 for profile in self.profiles:
801 if profile.nlayers >= nlayers:
802 r_container.append(profile)
804 return r_container
806 def selectMaxLayers(self, nlayers):
807 '''
808 Select profiles with more than ``nlayers``.
810 :param nlayers: Maximum number of layers
811 :type nlayers: int
812 :returns: Selected profiles
813 :rtype: :class:`CrustDB`
814 '''
815 r_container = self._emptyCopy()
816 logger.info('Selecting maximum %d layers...' % nlayers)
818 for profile in self.profiles:
819 if profile.nlayers <= nlayers:
820 r_container.append(profile)
822 return r_container
824 def selectMinDepth(self, depth):
825 '''
826 Select profiles describing layers deeper than ``depth``.
828 :param depth: Minumum depth in [m]
829 :type depth: float
830 :returns: Selected profiles
831 :rtype: :class:`CrustDB`
832 '''
833 r_container = self._emptyCopy()
834 logger.info('Selecting minimum depth %f m...' % depth)
836 for profile in self.profiles:
837 if profile.d.max() >= depth:
838 r_container.append(profile)
839 return r_container
841 def selectMaxDepth(self, depth):
842 '''
843 Select profiles describing layers shallower than ``depth``.
845 :param depth: Maximum depth in [m]
846 :type depth: float
847 :returns: Selected profiles
848 :rtype: :class:`CrustDB`
849 '''
850 r_container = self._emptyCopy()
851 logger.info('Selecting maximum depth %f m...' % depth)
853 for profile in self.profiles:
854 if profile.d.max() <= depth:
855 r_container.append(profile)
856 return r_container
858 def selectVp(self):
859 '''
860 Select profiles describing P wave velocity.
862 :returns Selected profiles
863 :rtype: :class:`CrustDB`
864 '''
865 r_container = self._emptyCopy()
866 logger.info('Selecting profiles providing Vp...')
868 for profile in self.profiles:
869 if not num.all(num.isnan(profile.vp)):
870 r_container.append(profile)
871 return r_container
873 def selectVs(self):
874 '''
875 Select profiles describing P wave velocity.
877 :returns: Selected profiles
878 :rtype: :class:`CrustDB`
879 '''
880 r_container = self._emptyCopy()
881 logger.info('Selecting profiles providing Vs...')
883 for profile in self.profiles:
884 if not num.all(num.isnan(profile.vs)):
885 r_container.append(profile)
886 return r_container
888 def _emptyCopy(self):
889 r_container = CrustDB(parent=self)
890 r_container.name = self.name
891 return r_container
893 def exportCSV(self, filename=None):
894 '''
895 Export as CSV file.
897 :param filename: Export filename
898 :type filename: str
899 '''
900 with open(filename, 'w') as file:
901 file.write('# uid, Lat, Lon, vp, vs, H, Depth, Reference\n')
902 for profile in self.profiles:
903 file.write(profile._csv())
905 def exportYAML(self, filename=None):
906 '''
907 Export as YAML file.
909 :param filename: Export filename
910 :type filename: str
911 '''
912 with open(filename, 'w') as file:
913 for profile in self.profiles:
914 file.write(profile.__str__())
916 @classmethod
917 def readDatabase(cls, database_file):
918 db = cls()
919 CrustDB._read(db, database_file)
920 return db
922 @staticmethod
923 def _getRepositoryDatabase():
924 from pyrocko import config
926 name = path.basename(db_url)
927 data_path = path.join(config.config().crustdb_dir, name)
928 if not path.exists(data_path):
929 from pyrocko import util
930 util.download_file(db_url, data_path, None, None)
932 return data_path
934 def _read(self, database_file):
935 '''
936 Reads in the the GSN database file and puts it in CrustDB.
938 File format:
940 uid lat/lon vp vs hc depth
941 2 29.76N 2.30 .00 2.00 .00 s 25.70 .10 .00 NAC-CO 5 U
942 96.31W 3.94 .00 5.30 2.00 s 33.00 MCz 39.00 61C.3 EXC
943 5.38 .00 12.50 7.30 c
944 6.92 .00 13.20 19.80 c
945 8.18 .00 .00 33.00 m
947 3 34.35N 3.00 .00 3.00 .00 s 35.00 1.60 .00 NAC-BR 4 R
948 117.83W 6.30 .00 16.50 3.00 38.00 MCz 55.00 63R.1 ORO
949 7.00 .00 18.50 19.50
950 7.80 .00 .00 38.00 m
953 :param database_file: path to database file, type string
955 '''
957 def get_empty_record():
958 meta = {
959 'uid': num.nan,
960 'geographical_location': None,
961 'geological_province': None,
962 'geological_age': None,
963 'elevation': num.nan,
964 'heatflow': num.nan,
965 'measurement_method': None,
966 'publication_reference': None
967 }
968 # vp, vs, h, d, lat, lon, meta
969 return (num.zeros(128, dtype=num.float32),
970 num.zeros(128, dtype=num.float32),
971 num.zeros(128, dtype=num.float32),
972 num.zeros(128, dtype=num.float32),
973 0., 0., meta)
975 nlayers = []
977 def add_record(vp, vs, h, d, lat, lon, meta, nlayer):
978 if nlayer == 0:
979 return
980 self.append(VelocityProfile(
981 vp=vp[:nlayer] * km,
982 vs=vs[:nlayer] * km,
983 h=h[:nlayer] * km,
984 d=d[:nlayer] * km,
985 lat=lat, lon=lon,
986 **meta))
987 nlayers.append(nlayer)
989 vp, vs, h, d, lat, lon, meta = get_empty_record()
990 ilayer = 0
991 with open(database_file, 'r') as database:
992 for line, dbline in enumerate(database):
993 if dbline.isspace():
994 if not len(d) == 0:
995 add_record(vp, vs, h, d, lat, lon, meta, ilayer)
996 if not len(vp) == len(h):
997 raise DatabaseError(
998 'Inconsistent database, check line %d!\n\tDebug: '
999 % line, lat, lon, vp, vs, h, d, meta)
1001 vp, vs, h, d, lat, lon, meta = get_empty_record()
1002 ilayer = 0
1003 else:
1004 try:
1005 if ilayer == 0:
1006 lat = float(dbline[8:13])
1007 if dbline[13] == b'S':
1008 lat = -lat
1009 # Additional meta data
1010 meta['uid'] = int(dbline[0:6])
1011 meta['elevation'] = float(dbline[52:57])
1012 meta['heatflow'] = float(dbline[58:64])
1013 if meta['heatflow'] == 0.:
1014 meta['heatflow'] = None
1015 meta['geographical_location'] =\
1016 dbline[66:72].strip()
1017 meta['measurement_method'] = dbline[77]
1018 if ilayer == 1:
1019 lon = float(dbline[7:13])
1020 if dbline[13] == b'W':
1021 lon = -lon
1022 # Additional meta data
1023 meta['geological_age'] = dbline[54:58].strip()
1024 meta['publication_reference'] =\
1025 dbline[66:72].strip()
1026 meta['geological_province'] = dbline[74:78].strip()
1027 try:
1028 vp[ilayer] = dbline[17:21]
1029 vs[ilayer] = dbline[23:27]
1030 h[ilayer] = dbline[28:34]
1031 d[ilayer] = dbline[35:41]
1032 except ValueError:
1033 pass
1034 except ValueError:
1035 logger.warning(
1036 'Could not interpret line %d, skipping\n%s' %
1037 (line, dbline))
1038 while not database.readlines():
1039 pass
1040 vp, vs, h, d, lat, lon, meta = get_empty_record()
1041 ilayer += 1
1042 # Append last profile
1043 add_record(vp, vs, h, d, lat, lon, meta, ilayer)
1044 logger.info('Loaded %d profiles from %s' %
1045 (self.nprofiles, database_file))