|
# http://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg---------- '''Interface to use CRUST2.0 model by Laske, Masters and Reif.
All functions defined in this module return SI units (m, m/s, kg/m^3).
.. note:: Please refer to the REM web site if you use this model:
http://igppweb.ucsd.edu/~gabi/rem.html
or
Bassin, C., Laske, G. and Masters, G., The Current Limits of Resolution for Surface Wave Tomography in North America, EOS Trans AGU, 81, F897, 2000. A description of CRUST 5.1 can be found in: Mooney, Laske and Masters, Crust 5.1: a global crustal model at 5x5 degrees, JGR, 103, 727-747, 1998.
Usage -----
::
>>> from pyrocko import crust2x2 >>> p = crust2x2.get_profile(10., 20.) >>> print p type, name: G2, Archean 0.5 km seds. elevation: 529 crustal thickness: 38500 average vp, vs, rho: 6460.7 3665.1 2867.5 mantle ave. vp, vs, rho: 8200 4700 3400
0 3810 1940 920 ice 0 1500 0 1020 water 500 2500 1200 2100 soft sed. 0 4000 2100 2400 hard sed. 12500 6200 3600 2800 upper crust 13000 6400 3600 2850 middle crust 13000 6800 3800 2950 lower crust >>> print p.get_weeded() [[ 0. 500. 500. 13000. 13000. 26000. 26000. 39000. 39000.] [ 2500. 2500. 6200. 6200. 6400. 6400. 6800. 6800. 8200.] [ 1200. 1200. 3600. 3600. 3600. 3600. 3800. 3800. 4700.] [ 2100. 2100. 2800. 2800. 2850. 2850. 2950. 2950. 3400.]]
Constants ---------
============== ============== Layer id Layer name ============== ============== LICE ice LWATER water LSOFTSED soft sediments LHARDSED hard sediments LUPPERCRUST upper crust LMIDDLECRUST middle crust LLOWERCRUST lower crust LBELOWCRUST below crust ============== ==============
Contents -------- ''' from __future__ import absolute_import, division
from builtins import range
import os import copy import math from io import StringIO
import numpy as num
LICE, LWATER, LSOFTSED, LHARDSED, LUPPERCRUST, LMIDDLECRUST, \ LLOWERCRUST, LBELOWCRUST = list(range(8))
class Crust2Profile(object): '''Representation of a CRUST2.0 key profile.'''
layer_names = ( 'ice', 'water', 'soft sed.', 'hard sed.', 'upper crust', 'middle crust', 'lower crust', 'mantle')
def __init__(self, ident, name, vp, vs, rho, thickness, elevation): self._ident = ident self._name = name self._vp = vp self._vs = vs self._rho = rho self._thickness = thickness self._crustal_thickness = None self._elevation = elevation
def get_weeded(self, include_waterlayer=False): '''Get layers used in the profile.
:param include_waterlayer: include water layer if ``True``. Default is ``False``
:returns: NumPy array with rows ``depth``, ``vp``, ``vs``, ``density`` ''' depth = 0. layers = [] for ilayer, thickness, vp, vs, rho in zip( range(8), self._thickness, self._vp[:-1], self._vs[:-1], self._rho[:-1]):
if thickness == 0.0: continue
if not include_waterlayer and ilayer == LWATER: continue
layers.append([depth, vp, vs, rho]) layers.append([depth+thickness, vp, vs, rho]) depth += thickness
layers.append([ depth, self._vp[LBELOWCRUST], self._vs[LBELOWCRUST], self._rho[LBELOWCRUST]])
return num.array(layers).T
def get_layer(self, ilayer): '''Get parameters for a layer.
:param ilayer: id of layer :returns: thickness, vp, vs, density '''
if ilayer == LBELOWCRUST: thickness = num.Inf else: thickness = self._thickness[ilayer]
return thickness, self._vp[ilayer], self._vs[ilayer], self._rho[ilayer]
def set_elevation(self, elevation): self._elevation = elevation
def set_layer_thickness(self, ilayer, thickness): self._thickness[ilayer] = thickness
def elevation(self): return self._elevation
def __str__(self):
vvp, vvs, vrho, vthi = self.averages()
return '''type, name: %s, %s elevation: %15.5g crustal thickness: %15.5g average vp, vs, rho: %15.5g %15.5g %15.5g mantle ave. vp, vs, rho: %15.5g %15.5g %15.5g
%s''' % (self._ident, self._name, self._elevation, vthi, vvp, vvs, vrho, self._vp[LBELOWCRUST], self._vs[LBELOWCRUST], self._rho[LBELOWCRUST], '\n'.join([ '%15.5g %15.5g %15.5g %15.5g %s' % x for x in zip( self._thickness, self._vp[:-1], self._vs[:-1], self._rho[:-1], Crust2Profile.layer_names)]))
def crustal_thickness(self): '''Get total crustal thickness
Takes into account ice layer. Does not take into account water layer. '''
return num.sum(self._thickness[3:]) + self._thickness[LICE]
def averages(self): '''Get crustal averages for vp, vs and density and total crustal thickness,
Takes into account ice layer. Does not take into account water layer. '''
vthi = self.crustal_thickness() vvp = num.sum(self._thickness[3:] / self._vp[3:-1]) + \ self._thickness[LICE] / self._vp[LICE] vvs = num.sum(self._thickness[3:] / self._vs[3:-1]) + \ self._thickness[LICE] / self._vs[LICE] vrho = num.sum(self._thickness[3:] * self._rho[3:-1]) + \ self._thickness[LICE] * self._rho[LICE]
vvp = vthi / vvp vvs = vthi / vvs vrho = vrho / vthi
return vvp, vvs, vrho, vthi
def _sa2arr(sa): return num.array([float(x) for x in sa], dtype=num.float)
def _wrap(x, mi, ma): if mi <= x and x <= ma: return x
return x - math.floor((x-mi)/(ma-mi)) * (ma-mi)
def _clip(x, mi, ma): return min(max(mi, x), ma)
class Crust2(object): '''Access CRUST2.0 model.
:param directory: Directory with the data files which contain the CRUST2.0 model data. If this is set to ``None``, builtin CRUST2.0 files are used. '''
fn_keys = 'CNtype2_key.txt' fn_elevation = 'CNelevatio2.txt' fn_map = 'CNtype2.txt'
nlo = 180 nla = 90
_instance = None
def __init__(self, directory=None):
self.profile_keys = [] self._directory = directory self._typemap = None self._load_crustal_model()
def get_profile(self, *args, **kwargs): ''' Get crustal profile at a specific location or raw profile for given key.
Get profile for location ``(lat, lon)``, or raw profile for given string key.
:rtype: instance of :py:class:`Crust2Profile` '''
lat = kwargs.pop('lat', None) lon = kwargs.pop('lon', None)
if len(args) == 2: lat, lon = args
if lat is not None and lon is not None: return self._typemap[self._indices(float(lat), float(lon))] else: return self._raw_profiles[args[0]]
def _indices(self, lat, lon): lat = _clip(lat, -90., 90.) lon = _wrap(lon, -180., 180.) dlo = 360./Crust2.nlo dla = 180./Crust2.nla cola = 90.-lat ilat = _clip(int(cola/dla), 0, Crust2.nla-1) ilon = int((lon+180.)/dlo) % Crust2.nlo return ilat, ilon
def _load_crustal_model(self):
if self._directory is not None: path_keys = os.path.join(self._directory, Crust2.fn_keys) f = open(path_keys, 'r') else: from .crust2x2_data import decode, type2_key, type2, elevation
f = StringIO(decode(type2_key))
# skip header for i in range(5): f.readline()
profiles = {} while True: line = f.readline() if not line: break ident, name = line.split(None, 1) line = f.readline() vp = _sa2arr(line.split()) * 1000. line = f.readline() vs = _sa2arr(line.split()) * 1000. line = f.readline() rho = _sa2arr(line.split()) * 1000. line = f.readline() toks = line.split() thickness = _sa2arr(toks[:-2]) * 1000.
assert ident not in profiles
profiles[ident] = Crust2Profile( ident.strip(), name.strip(), vp, vs, rho, thickness, 0.0)
f.close()
self._raw_profiles = profiles self.profile_keys = sorted(profiles.keys())
if self._directory is not None: path_map = os.path.join(self._directory, Crust2.fn_map) f = open(path_map, 'r') else: f = StringIO(decode(type2))
f.readline() # header
amap = {} for ila, line in enumerate(f): keys = line.split()[1:] for ilo, key in enumerate(keys): amap[ila, ilo] = copy.deepcopy(profiles[key])
f.close()
if self._directory is not None: path_elevation = os.path.join(self._directory, Crust2.fn_elevation) f = open(path_elevation, 'r')
else: f = StringIO(decode(elevation))
f.readline() for ila, line in enumerate(f): for ilo, s in enumerate(line.split()[1:]): p = amap[ila, ilo] p.set_elevation(float(s)) if p.elevation() < 0.: p.set_layer_thickness(LWATER, -p.elevation())
f.close()
self._typemap = amap
@staticmethod def instance(): '''Get the global default Crust2 instance.'''
if Crust2._instance is None: Crust2._instance = Crust2()
return Crust2._instance
def get_profile_keys(): '''Get list of all profile keys.'''
crust2 = Crust2.instance() return list(crust2.profile_keys)
def get_profile(*args, **kwargs): '''Get Crust2x2 profile for given location or profile key.
Get profile for (lat,lon) or raw profile for given string key. '''
crust2 = Crust2.instance() return crust2.get_profile(*args, **kwargs)
def plot_crustal_thickness(crust2=None, filename='crustal_thickness.pdf'): ''' Create a quick and dirty plot of the crustal thicknesses defined in CRUST2.0. '''
if crust2 is None: crust2 = Crust2.instance()
def func(lat, lon): return crust2.get_profile(lat, lon).crustal_thickness(),
plot(func, filename, zscaled_unit='km', zscaled_unit_factor=0.001)
def plot_vp_belowcrust(crust2=None, filename='vp_below_crust.pdf'): ''' Create a quick and dirty plot of vp below the crust, as defined in CRUST2.0. '''
if crust2 is None: crust2 = Crust2.instance()
def func(lat, lon): return crust2.get_profile(lat, lon).get_layer(LBELOWCRUST)[1]
plot(func, filename, zscaled_unit='km/s', zscaled_unit_factor=0.001)
def plot(func, filename, **kwargs): nlats, nlons = 91, 181 lats = num.linspace(-90., 90., nlats) lons = num.linspace(-180., 180., nlons)
vecfunc = num.vectorize(func, [num.float]) latss, lonss = num.meshgrid(lats, lons) thickness = vecfunc(latss, lonss)
from pyrocko.plot import gmtpy cm = gmtpy.cm marg = (1.5*cm, 2.5*cm, 1.5*cm, 1.5*cm) p = gmtpy.Simple( width=20*cm, height=10*cm, margins=marg, with_palette=True, **kwargs)
p.density_plot(gmtpy.tabledata(lons, lats, thickness.T)) p.save(filename) |